#sequencing data used to generat OTU tables via Qiime2 version 2018.11, can be found under NCBI SRA accession: PRJNA615043
getwd()
## [1] "/Users/martinylab/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020"
setwd("/Users/martinylab/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020")
library(BiodiversityR)
## Loading required package: tcltk
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
## Loading required package: vegan3d
## Loading required package: rgl
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
## BiodiversityR 2.11-3: Use command BiodiversityRGUI() to launch the Graphical User Interface;
## to see changes use BiodiversityRGUI(changeLog=TRUE, backward.compatibility.messages=TRUE)
library(biomformat)
library(car)
## Loading required package: carData
library(data.table)
library(doBy)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
##
## between, first, last
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(devtools)
## Loading required package: usethis
##
## Attaching package: 'devtools'
## The following object is masked from 'package:permute':
##
## check
library(EcolUtils)
library(emmeans)
##
## Attaching package: 'emmeans'
## The following object is masked from 'package:devtools':
##
## test
library(forcats)
library(gapminder)
library(ggplot2)
library(ggpubr)
library(lmerTest)
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(lsmeans)
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:dplyr':
##
## intersect, setdiff, union
## The following objects are masked from 'package:data.table':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(multcompView)
library(purrr)
##
## Attaching package: 'purrr'
## The following object is masked from 'package:data.table':
##
## transpose
## The following object is masked from 'package:car':
##
## some
library(qiime2R)
library(readxl)
library(rcompanion)
library(Rmisc)
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following object is masked from 'package:purrr':
##
## compact
## The following object is masked from 'package:ggpubr':
##
## mutate
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
library(stringr)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ tibble 3.0.1 ✓ readr 1.3.1
## ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x plyr::arrange() masks dplyr::arrange()
## x lubridate::as.difftime() masks base::as.difftime()
## x dplyr::between() masks data.table::between()
## x readr::col_factor() masks scales::col_factor()
## x plyr::compact() masks purrr::compact()
## x plyr::count() masks dplyr::count()
## x lubridate::date() masks base::date()
## x scales::discard() masks purrr::discard()
## x tidyr::expand() masks Matrix::expand()
## x plyr::failwith() masks dplyr::failwith()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks data.table::first()
## x lubridate::hour() masks data.table::hour()
## x plyr::id() masks dplyr::id()
## x lubridate::intersect() masks base::intersect()
## x lubridate::isoweek() masks data.table::isoweek()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks data.table::last()
## x lubridate::mday() masks data.table::mday()
## x lubridate::minute() masks data.table::minute()
## x lubridate::month() masks data.table::month()
## x plyr::mutate() masks ggpubr::mutate(), dplyr::mutate()
## x tidyr::pack() masks Matrix::pack()
## x lubridate::quarter() masks data.table::quarter()
## x dplyr::recode() masks car::recode()
## x plyr::rename() masks dplyr::rename()
## x lubridate::second() masks data.table::second()
## x lubridate::setdiff() masks base::setdiff()
## x purrr::some() masks car::some()
## x plyr::summarise() masks dplyr::summarise()
## x plyr::summarize() masks dplyr::summarize()
## x purrr::transpose() masks data.table::transpose()
## x lubridate::union() masks base::union()
## x tidyr::unpack() masks Matrix::unpack()
## x lubridate::wday() masks data.table::wday()
## x lubridate::week() masks data.table::week()
## x lubridate::yday() masks data.table::yday()
## x lubridate::year() masks data.table::year()
library(vegan)
library(xlsx)
library(yaml)
###Bacteria###
unzip(zipfile = "table_filtered.16s.qza")
hdf5_biom.16s <- read_hdf5_biom("/Users/martinylab/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/1a5aa860-8bb9-46e9-8a69-2499a6e1cc21/data/feature-table.biom")
write_biom(hdf5_biom.16s, "formatted_biom")
json_biom.16s <- read_biom("formatted_biom")
otu.table.16s <- as.data.frame(as.matrix(biom_data(json_biom.16s)))
write.table(otu.table.16s, file = "otu.table.filtered.16s.tsv", quote=FALSE, sep='\t', col.names = NA)
###Fungi###
hdf5_biom.ITS <- read_hdf5_biom("~/Google_drive/my_projects/lrgce/lrgce_seq_v2/ITS/forward_final/exported/feature-table.biom")
write_biom(hdf5_biom.ITS, "formatted_biom")
json_biom.ITS <- read_biom("formatted_biom")
otu.table.ITS <- as.data.frame(as.matrix(biom_data(json_biom.ITS)))
write.table(otu.table.ITS, file = "otu.table.filtered.ITS.tsv", quote=FALSE, sep='\t', col.names = NA)
###Plant###
plant.2015.fg <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/plant.2015.fg.xlsx")
otu.plant.2015 <- as.data.frame(plant.2015.fg, col.names = TRUE, row.names = TRUE)
otu.plant.2015 [1:5 , 1:6]
## sample_id Bare Litter AVspp BRDI BRHO
## 1 G05LRN_2015 7.843137 4.901961 0 0.0000000 0
## 2 G05RRX_2015 11.764706 12.745098 0 0.0000000 0
## 3 G07LRX_2015 4.901961 2.941176 0 0.0000000 0
## 4 G07RRN_2015 15.686275 5.882353 0 6.8627451 0
## 5 G14LRN_2015 13.725490 14.705882 0 0.9803922 0
otu.plant.2015.1 <- otu.plant.2015[-1]
row.names(otu.plant.2015.1) <- otu.plant.2015$sample_id
otu.plant.2015.1[1:5 , 1:4]
## Bare Litter AVspp BRDI
## G05LRN_2015 7.843137 4.901961 0 0.0000000
## G05RRX_2015 11.764706 12.745098 0 0.0000000
## G07LRX_2015 4.901961 2.941176 0 0.0000000
## G07RRN_2015 15.686275 5.882353 0 6.8627451
## G14LRN_2015 13.725490 14.705882 0 0.9803922
combo.16s <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/combo.16s.xlsx")
metadata.16s <- as.data.frame(combo.16s, col.names = TRUE, row.names = TRUE)
metadata.16s[1:5 , 1:6]
## Sample_id vegtype_plot_treatments Ecosystem Plot Block Precipitation
## 1 14RRXAugust201616S G14RRX Grass 14 6 Reduced
## 2 18LXNAugust201616S G18LXN Grass 18 6 Ambient
## 3 18RXXAugust201616S G18RXX Grass 18 6 Ambient
## 4 20LRXAugust201616S G20LRX Grass 20 8 Reduced
## 5 22LXXAugust201616S G22LXX Grass 22 8 Ambient
metadata.16s.1 <- metadata.16s[-1]
row.names(metadata.16s.1) <- metadata.16s$Sample_id
sapply(metadata.16s.1, class)
## vegtype_plot_treatments Ecosystem
## "character" "character"
## Plot Block
## "numeric" "numeric"
## Precipitation Nitrogen
## "character" "character"
## eco.nit eco.prec
## "character" "character"
## Collection_date season
## "character" "character"
## eco.prec.seas Rainfall(in)
## "character" "numeric"
## Rainfall(cm) Rainfall_corrected_by_treatment(cm)
## "numeric" "numeric"
## shannon richness
## "numeric" "numeric"
## pielou_e average_sobs
## "numeric" "numeric"
metadata.16s.1[1:5 , 1:6]
## vegtype_plot_treatments Ecosystem Plot Block Precipitation
## 14RRXAugust201616S G14RRX Grass 14 6 Reduced
## 18LXNAugust201616S G18LXN Grass 18 6 Ambient
## 18RXXAugust201616S G18RXX Grass 18 6 Ambient
## 20LRXAugust201616S G20LRX Grass 20 8 Reduced
## 22LXXAugust201616S G22LXX Grass 22 8 Ambient
## Nitrogen
## 14RRXAugust201616S Ambient
## 18LXNAugust201616S Added
## 18RXXAugust201616S Ambient
## 20LRXAugust201616S Ambient
## 22LXXAugust201616S Ambient
combo.ITS <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/combo.ITS.xlsx")
metadata.ITS <- as.data.frame(combo.ITS, col.names = TRUE, row.names = TRUE)
metadata.ITS[1:5 , 1:6]
## sample vegtype_plot_treatments Plot Block Ecosystem Precipitation
## 1 14RRXAugust2016ITS G14RRX 14 6 Grass Reduced
## 2 18LXNAugust2016ITS G18LXN 18 6 Grass Ambient
## 3 18RXXAugust2016ITS G18RXX 18 6 Grass Ambient
## 4 20LRXAugust2016ITS G20LRX 20 8 Grass Reduced
## 5 22LXXAugust2016ITS G22LXX 22 8 Grass Ambient
metadata.ITS.1 <- metadata.ITS[-1]
row.names(metadata.ITS.1) <- metadata.ITS$sample
sapply(metadata.ITS.1, class)
## vegtype_plot_treatments Plot
## "character" "numeric"
## Block Ecosystem
## "numeric" "character"
## Precipitation Nitrogen
## "character" "character"
## eco.nit eco.prec
## "character" "character"
## Collection_date season
## "character" "character"
## Rainfall(in) Rainfall(cm)
## "numeric" "numeric"
## Rainfall_corrected_by_treatment(cm) shannon
## "numeric" "numeric"
## richness pielou_e
## "numeric" "numeric"
## average_sobs
## "numeric"
metadata.ITS.1[1:5 , 1:6]
## vegtype_plot_treatments Plot Block Ecosystem Precipitation
## 14RRXAugust2016ITS G14RRX 14 6 Grass Reduced
## 18LXNAugust2016ITS G18LXN 18 6 Grass Ambient
## 18RXXAugust2016ITS G18RXX 18 6 Grass Ambient
## 20LRXAugust2016ITS G20LRX 20 8 Grass Reduced
## 22LXXAugust2016ITS G22LXX 22 8 Grass Ambient
## Nitrogen
## 14RRXAugust2016ITS Ambient
## 18LXNAugust2016ITS Added
## 18RXXAugust2016ITS Ambient
## 20LRXAugust2016ITS Ambient
## 22LXXAugust2016ITS Ambient
plant.combo.meta.2015 <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/plant.combo.meta_2015.xlsx")
metadata.plant.2015 <- as.data.frame(plant.combo.meta.2015, col.names = TRUE, row.names = TRUE)
metadata.plant.2015[1:5 , 1:6]
## sample_id plot block year ecosystem nitrogen
## 1 G05LRN_2015 G05LRN 2 2015 grass added
## 2 G05RRX_2015 G05RRX 2 2015 grass ambient
## 3 G07LRX_2015 G07LRX 3 2015 grass ambient
## 4 G07RRN_2015 G07RRN 3 2015 grass added
## 5 G14LRN_2015 G14LRN 6 2015 grass added
metadata.plant.2015.1 <- metadata.plant.2015[-1]
row.names(metadata.plant.2015.1) <- metadata.plant.2015$sample_id
sapply(metadata.plant.2015.1, class)
## plot block year ecosystem
## "character" "numeric" "numeric" "character"
## nitrogen precipitation rainfall(cm) eco.prec
## "character" "character" "numeric" "character"
## eco.nit NonNativeForbCover NonNativeGrassCover NativeForbCover
## "character" "numeric" "numeric" "numeric"
## NativeGrassCover NativeShrubCover Bare Litter
## "numeric" "numeric" "numeric" "numeric"
## total_grass_cover TotalNonNativeCover TotalNativeCover
## "numeric" "numeric" "numeric"
metadata.plant.2015.1[1:5 , 1:6]
## plot block year ecosystem nitrogen precipitation
## G05LRN_2015 G05LRN 2 2015 grass added reduced
## G05RRX_2015 G05RRX 2 2015 grass ambient reduced
## G07LRX_2015 G07LRX 3 2015 grass ambient reduced
## G07RRN_2015 G07RRN 3 2015 grass added reduced
## G14LRN_2015 G14LRN 6 2015 grass added reduced
taxonomy.16s <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/taxonomy_16s.xlsx")
taxaid.16s <- as.data.frame(taxonomy.16s, col.names = TRUE, row.names = TRUE)
taxaid.16s[1:2 , 1:2]
## OTUID
## 1 9ccef5c0259d44ee21a37f1838ef585b
## 2 6f301efd81a9f8bfb9c3d39ff3700c8a
## taxonomy
## 1 k__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Actinomycetales; f__Streptomycetaceae; g__Streptomyces; s__reticuliscabiei
## 2 k__Bacteria; p__Bacteroidetes; c__Sphingobacteriia; o__Sphingobacteriales; f__; g__; s__
taxaid.16s.1 <-taxaid.16s[-1]
row.names(taxaid.16s.1) <- taxaid.16s$OTUID
sapply(taxaid.16s.1, class)
## taxonomy confidence
## "character" "numeric"
taxaid.16s.1[1:2 , 1:2]
## taxonomy
## 9ccef5c0259d44ee21a37f1838ef585b k__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Actinomycetales; f__Streptomycetaceae; g__Streptomyces; s__reticuliscabiei
## 6f301efd81a9f8bfb9c3d39ff3700c8a k__Bacteria; p__Bacteroidetes; c__Sphingobacteriia; o__Sphingobacteriales; f__; g__; s__
## confidence
## 9ccef5c0259d44ee21a37f1838ef585b 0.9828906
## 6f301efd81a9f8bfb9c3d39ff3700c8a 0.9999996
taxonomy.ITS <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/taxonomy_ITS.xlsx")
taxaid.ITS <- as.data.frame(taxonomy.ITS, col.names = TRUE, row.names = TRUE)
taxaid.ITS[1:2 , 1:2]
## otu_id
## 1 7e80c7e5f5813d51bcb914eae7555273
## 2 1094975b0bdad9f93c92f4b3b2e6df15
## taxon
## 1 k__Fungi;p__Ascomycota
## 2 k__Fungi;p__Ascomycota;c__Leotiomycetes;o__Helotiales;f__Helotiaceae;g__Articulospora;s__Articulospora_proliferata
taxaid.ITS.1 <-taxaid.ITS[-1]
row.names(taxaid.ITS.1) <- taxaid.ITS$otu_id
sapply(taxaid.ITS.1, class)
## taxon confidence
## "character" "numeric"
taxaid.ITS.1[1:2 , 1:2]
## taxon
## 7e80c7e5f5813d51bcb914eae7555273 k__Fungi;p__Ascomycota
## 1094975b0bdad9f93c92f4b3b2e6df15 k__Fungi;p__Ascomycota;c__Leotiomycetes;o__Helotiales;f__Helotiaceae;g__Articulospora;s__Articulospora_proliferata
## confidence
## 7e80c7e5f5813d51bcb914eae7555273 0.8120831
## 1094975b0bdad9f93c92f4b3b2e6df15 0.9517936
###Bacteria###
otu.plus.taxaid.16s <- as.data.frame(merge(taxaid.16s.1, otu.table.16s, by.x = "row.names", by.y = "row.names"))
otu.plus.taxaid.filtered.16s <- otu.plus.taxaid.16s[!grepl("chloroplast", otu.plus.taxaid.16s$`taxon`),]
otu.plus.taxaid.filtered.16s <- otu.plus.taxaid.filtered.16s[!grepl("mitochondria",otu.plus.taxaid.filtered.16s$`taxon`),]
write.table(otu.plus.taxaid.filtered.16s, file = "otu.taxaid.filtered.16s.tsv", quote=FALSE, sep='\t', col.names = NA)
row.names(otu.plus.taxaid.filtered.16s) <- otu.plus.taxaid.filtered.16s$Row.names
otu.clean.16s <- as.data.frame(t(otu.plus.taxaid.filtered.16s[,4:ncol(otu.plus.taxaid.filtered.16s)]))
write.table(otu.clean.16s, file = "otu.clean.16s.tsv", quote=FALSE, sep='\t', col.names = NA)
read.depth.barplot.16s <- barplot(sort(rowSums(otu.clean.16s)), ylim = c(0, max(rowSums(otu.clean.16s))),
xlim = c(0, NROW(otu.clean.16s)), col = "Blue", ylab = "Read Depth", xlab = "Sample")

unrarfied.curves.16s <-rarecurve(otu.clean.16s, step = 1000, label = FALSE)

set.seed(seed = 999)
rared.otu.16s <- as.data.frame((rrarefy.perm(otu.clean.16s, sample = 1090, n = 300, round.out = T)))#because rrarefy.perm randomly picks OTUs at a specified sequencing depth, the total number of OTUs may change with each run of the function.
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 2 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 3 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 4 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 5 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 6 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 7 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 8 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 9 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 10 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 11 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 12 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 13 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 14 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 15 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 16 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 17 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 18 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 19 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 20 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 21 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 22 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 23 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 24 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 25 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 26 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 27 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 28 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 29 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 30 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 31 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 32 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 33 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 34 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 35 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 36 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 37 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 38 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 39 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 40 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 41 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 42 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 43 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 44 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 45 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 46 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 47 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 48 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 49 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 50 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 51 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 52 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 53 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 54 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 55 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 56 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 57 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 58 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 59 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 60 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 61 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 62 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 63 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 64 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 65 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 66 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 67 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 68 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 69 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 70 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 71 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 72 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 73 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 74 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 75 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 76 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 77 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 78 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 79 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 80 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 81 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 82 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 83 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 84 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 85 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 86 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 87 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 88 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 89 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 90 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 91 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 92 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 93 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 94 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 95 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 96 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 97 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 98 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 99 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 100 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 101 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 102 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 103 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 104 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 105 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 106 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 107 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 108 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 109 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 110 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 111 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 112 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 113 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 114 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 115 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 116 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 117 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 118 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 119 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 120 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 121 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 122 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 123 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 124 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 125 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 126 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 127 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 128 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 129 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 130 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 131 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 132 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 133 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 134 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 135 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 136 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 137 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 138 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 139 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 140 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 141 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 142 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 143 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 144 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 145 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 146 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 147 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 148 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 149 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 150 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 151 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 152 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 153 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 154 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 155 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 156 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 157 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 158 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 159 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 160 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 161 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 162 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 163 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 164 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 165 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 166 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 167 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 168 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 169 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 170 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 171 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 172 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 173 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 174 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 175 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 176 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 177 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 178 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 179 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 180 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 181 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 182 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 183 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 184 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 185 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 186 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 187 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 188 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 189 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 190 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 191 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 192 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 193 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 194 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 195 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 196 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 197 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 198 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 199 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 200 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 201 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 202 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 203 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 204 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 205 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 206 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 207 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 208 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 209 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 210 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 211 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 212 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 213 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 214 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 215 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 216 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 217 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 218 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 219 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 220 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 221 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 222 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 223 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 224 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 225 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 226 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 227 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 228 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 229 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 230 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 231 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 232 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 233 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 234 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 235 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 236 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 237 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 238 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 239 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 240 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 241 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 242 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 243 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 244 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 245 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 246 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 247 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 248 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 249 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 250 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 251 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 252 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 253 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 254 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 255 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 256 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 257 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 258 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 259 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 260 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 261 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 262 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 263 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 264 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 265 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 266 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 267 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 268 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 269 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 270 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 271 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 272 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 273 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 274 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 275 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 276 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 277 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 278 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 279 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 280 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 281 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 282 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 283 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 284 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 285 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 286 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 287 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 288 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 289 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 290 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 291 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 292 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 293 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 294 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 295 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 296 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 297 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 298 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 299 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 300 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
warnings(rared.otu.16s)
rared.otu.16s <- as.data.frame(rared.otu.16s[rowSums(rared.otu.16s) >= 1090-(1090*.1), colSums(rared.otu.16s) >= 1])
rarefied.curves.16s <- rarecurve(rared.otu.16s, step = 300)

rared.otu.16s.t <- t(rared.otu.16s)
rared.otu.taxaid.16s <- merge(rared.otu.16s.t, taxaid.16s.1, by.x = "row.names", by.y = "row.names")
write.table(rared.otu.taxaid.16s, file = "rared.otu.taxaid.16s.tsv", quote=FALSE, sep='\t', col.names = NA)
###Fungi###
otu.plus.taxaid.ITS <- as.data.frame(merge(taxaid.ITS.1, otu.table.ITS, by.x = "row.names", by.y = "row.names"))
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.ITS[!grepl("k__Chromista", otu.plus.taxaid.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Rhizaria", otu.plus.taxaid.filtered.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Plantae", otu.plus.taxaid.filtered.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Alveolata", otu.plus.taxaid.filtered.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Metazoa", otu.plus.taxaid.filtered.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Fungi;p__unidentified", otu.plus.taxaid.filtered.ITS$`taxon`),]
write.table(otu.plus.taxaid.filtered.ITS, file = "otu.taxaid.filtered.ITS.tsv", quote=FALSE, sep='\t', col.names = NA)
row.names(otu.plus.taxaid.filtered.ITS) <- otu.plus.taxaid.filtered.ITS$Row.names
otu.clean.ITS <- as.data.frame(t(otu.plus.taxaid.filtered.ITS[,4:ncol(otu.plus.taxaid.filtered.ITS)]))
write.table(otu.clean.ITS, file = "otu.clean.ITS.tsv", quote=FALSE, sep='\t', col.names = NA)
read.depth.barplot.ITS <- barplot(sort(rowSums(otu.clean.ITS)), ylim = c(0, max(rowSums(otu.clean.ITS))),
xlim = c(0, NROW(otu.clean.ITS)), col = "Blue", ylab = "Read Depth", xlab = "Sample")

unrarfied.curves.ITS <-rarecurve(otu.clean.ITS, step = 1000, label = FALSE)

set.seed(seed = 999)
rared.otu.ITS <- as.data.frame((rrarefy.perm(otu.clean.ITS, sample = 1064, n = 300, round.out = T)))
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 2 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 3 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 4 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 5 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 6 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 7 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 8 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 9 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 10 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 11 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 12 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 13 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 14 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 15 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 16 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 17 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 18 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 19 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 20 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 21 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 22 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 23 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 24 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 25 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 26 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 27 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 28 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 29 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 30 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 31 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 32 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 33 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 34 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 35 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 36 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 37 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 38 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 39 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 40 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 41 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 42 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 43 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 44 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 45 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 46 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 47 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 48 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 49 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 50 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 51 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 52 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 53 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 54 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 55 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 56 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 57 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 58 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 59 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 60 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 61 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 62 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 63 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 64 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 65 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 66 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 67 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 68 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 69 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 70 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 71 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 72 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 73 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 74 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 75 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 76 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 77 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 78 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 79 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 80 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 81 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 82 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 83 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 84 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 85 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 86 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 87 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 88 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 89 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 90 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 91 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 92 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 93 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 94 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 95 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 96 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 97 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 98 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 99 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 100 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 101 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 102 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 103 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 104 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 105 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 106 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 107 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 108 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 109 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 110 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 111 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 112 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 113 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 114 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 115 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 116 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 117 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 118 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 119 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 120 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 121 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 122 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 123 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 124 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 125 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 126 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 127 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 128 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 129 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 130 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 131 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 132 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 133 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 134 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 135 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 136 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 137 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 138 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 139 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 140 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 141 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 142 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 143 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 144 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 145 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 146 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 147 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 148 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 149 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 150 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 151 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 152 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 153 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 154 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 155 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 156 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 157 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 158 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 159 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 160 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 161 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 162 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 163 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 164 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 165 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 166 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 167 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 168 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 169 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 170 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 171 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 172 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 173 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 174 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 175 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 176 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 177 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 178 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 179 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 180 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 181 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 182 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 183 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 184 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 185 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 186 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 187 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 188 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 189 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 190 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 191 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 192 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 193 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 194 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 195 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 196 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 197 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 198 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 199 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 200 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 201 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 202 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 203 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 204 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 205 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 206 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 207 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 208 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 209 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 210 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 211 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 212 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 213 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 214 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 215 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 216 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 217 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 218 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 219 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 220 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 221 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 222 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 223 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 224 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 225 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 226 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 227 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 228 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 229 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 230 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 231 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 232 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 233 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 234 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 235 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 236 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 237 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 238 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 239 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 240 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 241 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 242 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 243 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 244 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 245 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 246 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 247 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 248 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 249 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 250 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 251 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 252 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 253 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 254 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 255 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 256 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 257 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 258 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 259 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 260 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 261 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 262 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 263 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 264 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 265 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 266 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 267 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 268 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 269 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 270 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 271 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 272 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 273 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 274 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 275 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 276 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 277 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 278 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 279 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 280 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 281 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 282 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 283 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 284 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 285 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 286 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 287 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 288 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 289 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 290 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 291 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 292 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 293 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 294 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 295 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 296 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 297 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 298 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 299 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation 300 out of 300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
warnings(rared.otu.ITS)
rared.otu.ITS <- as.data.frame(rared.otu.ITS[rowSums(rared.otu.ITS) >= 1064-(1064*.1), colSums(rared.otu.ITS) >= 1])
rarefied.curves.ITS <- rarecurve(rared.otu.ITS, step = 300)

rared.otu.ITS.t <- t(rared.otu.ITS)
rared.otu.taxaid.ITS <- merge(rared.otu.ITS.t, taxaid.ITS.1, by.x = "row.names", by.y = "row.names")
write.table(rared.otu.taxaid.ITS, file = "rared.otu.taxaid.ITS.tsv", quote=FALSE, sep='\t', col.names = NA)
#Create function to read qza file into R
read_qza <- function(file, tmp, rm) {
if(missing(tmp)){tmp <- tempdir()}
if(missing(file)){stop("Path to artifact (.qza) not provided")}
if(missing(rm)){rm=TRUE} #remove the decompressed object from tmp
unzip(file, exdir=tmp)
unpacked<-unzip(file, exdir=tmp, list=TRUE)
artifact<-read_yaml(paste0(tmp,"/", paste0(gsub("/..+","", unpacked$Name[1]),"/metadata.yaml"))) #start by loading in the metadata not assuming it will be first file listed
artifact$contents<-data.frame(files=unpacked)
artifact$contents$size=sapply(paste0(tmp, "/", artifact$contents$files), file.size)
artifact$version=read.table(paste0(tmp,"/",artifact$uuid, "/VERSION"))
#if(sum(artifact$version$V2==c("2","4","2018.4.0"))!=3){warning("Artifact was not generated with Qiime2 2018.4, if data is not successfully imported, please report here github.com/jbisanz/qiime2R/issues")}#check version and throw warning if new format
#get data dependent on format
if(grepl("BIOMV", artifact$format)){
suppressWarnings(artifact$data<-as(biom_data(read_biom(paste0(tmp, "/", artifact$uui,"/data/feature-table.biom"))),"matrix")) #suppressing warning about \n characters
} else if (artifact$format=="NewickDirectoryFormat"){
artifact$data<-read.tree(paste0(tmp,"/",artifact$uuid,"/data/tree.nwk"))
} else if (artifact$format=="DistanceMatrixDirectoryFormat") {
artifact$data<-as.dist(read.table(paste0(tmp,"/", artifact$uuid, "/data/distance-matrix.tsv"), header=TRUE, row.names=1))
} else if (grepl("StatsDirFmt", artifact$format)) {
if(paste0(artifact$uuid, "/data/stats.csv") %in% artifact$contents$files.Name){artifact$data<-read.csv(paste0(tmp,"/", artifact$uuid, "/data/stats.csv"), header=TRUE, row.names=1)}
if(paste0(artifact$uuid, "/data/stats.tsv") %in% artifact$contents$files.Name){artifact$data<-read.table(paste0(tmp,"/", artifact$uuid, "/data/stats.tsv"), header=TRUE, row.names=1, sep='\t')} #can be tsv or csv
} else if (artifact$format=="TSVTaxonomyDirectoryFormat"){
artifact$data<-read.table(paste0(tmp,"/", artifact$uuid, "/data/taxonomy.tsv"), sep='\t', header=TRUE, quote="")
} else if (artifact$format=="OrdinationDirectoryFormat"){
linesplit<-suppressWarnings(readLines(paste0(tmp,"/", artifact$uuid, "/data/ordination.txt")))
linesplit<-linesplit[sapply(linesplit, function(x) x!="")]
for (i in 1:length(linesplit)){
if(grepl("^Eigvals\\t|^Proportion explained\\t|^Species\\t|^Site\\t|^Biplot\\t|^Site constraints\\t", linesplit[i])){
curfile=strsplit(linesplit[i],"\t")[[1]][1]
} else {
write(linesplit[i], paste0(tmp,"/", artifact$uuid, "/data/",curfile,".tmp"), append=TRUE)
}
}
for (outs in list.files(paste0(tmp,"/", artifact$uuid,"/data"), full.names = TRUE, pattern = "\\.tmp")){
NewLab<-gsub(" ", "", toTitleCase(gsub("\\.tmp", "", basename(outs))))
artifact$data[[NewLab]]<-read.table(outs,sep='\t', header=FALSE)
if(NewLab %in% c("Eigvals","ProportionExplained")){colnames(artifact$data[[NewLab]])<-paste0("PC",1:ncol(artifact$data[[NewLab]]))}
if(NewLab %in% c("Site","SiteConstraints")){colnames(artifact$data[[NewLab]])<-c("SampleID", paste0("PC",1:(ncol(artifact$data[[NewLab]])-1)))}
if(NewLab %in% c("Species")){colnames(artifact$data[[NewLab]])<-c("FeatureID", paste0("PC",1:(ncol(artifact$data[[NewLab]])-1)))}
}
artifact$data$Vectors<-artifact$data$Site #Rename Site to Vectors so this matches up with the syntax used in the tutorials
artifact$data$Site<-NULL
} else if (artifact$format=="DNASequencesDirectoryFormat") {
artifact$data<-readDNAStringSet(paste0(tmp,"/",artifact$uuid,"/data/dna-sequences.fasta"))
} else if (artifact$format=="AlignedDNASequencesDirectoryFormat") {
artifact$data<-readDNAMultipleAlignment(paste0(tmp,"/",artifact$uuid,"/data/aligned-dna-sequences.fasta"))
} else if (grepl("EMPPairedEndDirFmt|EMPSingleEndDirFmt|FastqGzFormat|MultiplexedPairedEndBarcodeInSequenceDirFmt|MultiplexedSingleEndBarcodeInSequenceDirFmt|PairedDNASequencesDirectoryFormat|SingleLanePerSamplePairedEndFastqDirFmt|SingleLanePerSampleSingleEndFastqDirFmt", artifact$format)) {
artifact$data<-data.frame(files=list.files(paste0(tmp,"/", artifact$uuid,"/data")))
artifact$data$size<-format(sapply(artifact$data$files, function(x){file.size(paste0(tmp,"/",artifact$uuid,"/data/",x))}, simplify = TRUE))
} else if (artifact$format=="AlphaDiversityDirectoryFormat") {
artifact$data<-read.table(paste0(tmp, "/", artifact$uuid, "/data/alpha-diversity.tsv"))
} else {
message("Format not supported, only a list of internal files and provenance is being imported.")
artifact$data<-list.files(paste0(tmp,"/",artifact$uuid, "/data"))
}
pfiles<-paste0(tmp,"/", grep("..+provenance/..+action.yaml", unpacked$Name, value=TRUE))
artifact$provenance<-lapply(pfiles, read_yaml)
names(artifact$provenance)<-grep("..+provenance/..+action.yaml", unpacked$Name, value=TRUE)
if(rm==TRUE){unlink(paste0(tmp,"/", artifact$uuid), recursive=TRUE)}
return(artifact)
}
###Bacteria###
taxonomy16s<-read_qza("taxonomy16s.qza", "./tmp", rm=T)
taxonomy16s$uuid
## [1] "67172321-d90d-436c-8328-b24d8743c424"
taxtable16s<-taxonomy16s$data %>% as.tibble() %>% separate(Taxon, sep="; ", c("Kingdom","Phylum","Class","Order","Family","Genus","Species")) #convert the table into a tabular split version
## Warning: `as.tibble()` is deprecated as of tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: Expected 7 pieces. Missing pieces filled with `NA` in 749 rows [3, 10,
## 12, 15, 16, 19, 24, 29, 34, 37, 40, 42, 45, 53, 55, 57, 60, 61, 62, 63, ...].
taxtable16s
## # A tibble: 1,866 x 9
## Feature.ID Kingdom Phylum Class Order Family Genus Species Confidence
## <fct> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 9ccef5c0259… k__Bact… p__Act… c__Ac… o__Ac… f__Str… g__St… s__ret… 0.983
## 2 6f301efd81a… k__Bact… p__Bac… c__Sp… o__Sp… f__ g__ s__ 1.00
## 3 9b1857d7d74… k__Bact… <NA> <NA> <NA> <NA> <NA> <NA> 0.960
## 4 5888bd387cc… k__Bact… p__Pro… c__De… o__My… f__ g__ s__ 0.798
## 5 f26c82353ae… k__Bact… p__Aci… c__[C… o__RB… f__Ell… g__ s__ 1.00
## 6 8e9ec770809… k__Bact… p__Bac… c__Fl… o__Fl… f__[We… g__Ch… s__ 1.00
## 7 6419cb3e5b3… k__Bact… p__Act… c__Ac… o__Ac… f__San… g__Sa… s__ 0.993
## 8 fb33e463293… k__Bact… p__Pro… c__Al… o__Sp… f__Sph… g__Sp… s__ 0.907
## 9 cadb3361e0a… k__Bact… p__Act… c__Ac… o__Ac… f__Kin… g__Ki… s__ 0.998
## 10 aa97c16a0a1… k__Bact… p__Pro… c__Be… o__Bu… f__Com… g__Li… <NA> 0.703
## # … with 1,856 more rows
taxtable16s <- as.data.frame(taxtable16s)
#remove p__ as a precursor to phylum name
taxtable16s$Phylum <- gsub("p__", "", taxtable16s$Phylum)
taxtable16s$Class <- gsub("c__", "", taxtable16s$Class)
taxtable16s$Order <- gsub("o__", "", taxtable16s$Order)
taxtable16s$Family <- gsub("f__", "", taxtable16s$Family)
taxtable16s$Genus <- gsub("g__", "", taxtable16s$Genus)
taxtable16s$Species <- gsub("s__", "", taxtable16s$Species)
#this will turn all empty cell into NA
taxtable16s <- taxtable16s %>%
mutate_at(vars(colnames(.)),
.funs = funs(ifelse(.=="", NA, as.character(.))))
## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.
#taxtable was imported using qiime2R, with already separated levels
row.names(taxtable16s) <- taxtable16s$Feature.ID #turn feature IDs into row.names
taxtable16s$Feature.ID <- NULL #remove column "Feature ID"
#merge OTU table with taxonomy table
otu.raref.taxonomy.16s <- as.data.frame(merge(taxtable16s, rared.otu.16s.t, by.x = "row.names", by.y = "row.names"))
row.names(otu.raref.taxonomy.16s) <- otu.raref.taxonomy.16s$Row.names #move column "Row.names" (first column) to row.names this was created during the merge
otu.raref.taxonomy.16s$Row.names <- NULL #remove column "Row.names"
otu.raref.taxonomy.16s$Confidence <- NULL # remove column "Confidence"
otu.raref.taxonomy.16s.ss <- subset(otu.raref.taxonomy.16s, select=c(8:195,1:7)) # move Taxa information to the end of all columns
write.table(otu.raref.taxonomy.16s.ss, "otu_raref_taxonomy.16s.csv", quote=FALSE, sep=',') #use this table to calculate top number of ESVs and their sequence contributions.
colSums(otu.raref.taxonomy.16s.ss[, c(1:188)]) # check if samples are RAREFIED!!
## 14LRNDecember201616S 14LRNJune201716S 14LRNMarch201716S
## 1090 1090 1097
## 14LRNMarch201816S 14RRXAugust201616S 14RRXDecember201616S
## 1090 1088 1088
## 14RRXDecember201716S 14RRXJune201716S 14RRXMarch201716S
## 1091 1090 1092
## 14RRXMarch201816S 14RRXSeptember201716S 18LXNAugust201616S
## 1093 1089 1087
## 18LXNDecember201616S 18LXNDecember201716S 18LXNJune201716S
## 1093 1089 1090
## 18LXNMarch201716S 18LXNMarch201816S 18LXNSeptember201716S
## 1091 1089 1091
## 18RXXAugust201616S 18RXXDecember201716S 18RXXJune201716S
## 1090 1091 1089
## 18RXXMarch201716S 18RXXMarch201816S 18RXXSeptember201716S
## 1089 1087 1090
## 20LRXAugust201616S 20LRXDecember201616S 20LRXDecember201716S
## 1090 1087 1090
## 20LRXJune201716S 20LRXMarch201716S 20LRXMarch201816S
## 1088 1091 1090
## 20LRXSeptember201716S 20RRNDecember201616S 20RRNDecember201716S
## 1089 1087 1091
## 20RRNMarch201716S 20RRNMarch201816S 22LXXAugust201616S
## 1092 1088 1091
## 22LXXDecember201616S 22LXXDecember201716S 22LXXMarch201716S
## 1092 1090 1092
## 22LXXMarch201816S 22RXNAugust201616S 22RXNDecember201616S
## 1089 1088 1091
## 22RXNDecember201716S 22RXNJune201716S 22RXNMarch201716S
## 1090 1088 1089
## 22RXNMarch201816S 22RXNSeptember201716S 25LRXAugust201616S
## 1088 1088 1089
## 25LRXDecember201616S 25LRXDecember201716S 25LRXMarch201716S
## 1086 1089 1090
## 25LRXMarch201816S 25LRXSeptember201716S 25RRNAugust201616S
## 1090 1092 1082
## 25RRNDecember201616S 25RRNDecember201716S 25RRNJune201716S
## 1089 1088 1091
## 25RRNMarch201716S 25RRNMarch201816S 25RRNSeptember201716S
## 1092 1092 1092
## 27LXNAugust201616S 27LXNDecember201616S 27LXNDecember201716S
## 1097 1096 1088
## 27LXNJune201716S 27LXNMarch201716S 27LXNMarch201816S
## 1091 1089 1091
## 27LXNSeptember201716S 27RXXAugust201616S 27RXXDecember201616S
## 1090 1086 1090
## 27RXXDecember201716S 27RXXJune201716S 27RXXMarch201716S
## 1087 1091 1087
## 27RXXMarch201816S 27RXXSeptember201716S 32LXNAugust201616S
## 1087 1090 1090
## 32LXNDecember201616S 32LXNMarch201716S 32LXNMarch201816S
## 1089 1088 1090
## 32LXNSeptember201716S 32RXXAugust201616S 32RXXDecember201616S
## 1097 1091 1090
## 32RXXDecember201716S 32RXXJune201716S 32RXXMarch201716S
## 1090 1092 1092
## 32RXXMarch201816S 32RXXSeptember201716S 35LRXAugust201616S
## 1089 1090 1089
## 35LRXDecember201616S 35LRXDecember201716S 35LRXJune201716S
## 1093 1091 1090
## 35LRXMarch201716S 35LRXMarch201816S 35LRXSeptember201716S
## 1090 1088 1090
## 35RRNAugust201616S 35RRNDecember201616S 35RRNMarch201716S
## 1090 1091 1090
## 35RRNMarch201816S 35RRNSeptember201716S 45LRXAugust201616S
## 1093 1091 1089
## 45LRXDecember201616S 45LRXJune201716S 45LRXMarch201716S
## 1089 1091 1091
## 45LRXMarch201816S 45LRXSeptember201716S 45RRNAugust201616S
## 1090 1089 1092
## 45RRNDecember201616S 45RRNDecember201716S 45RRNJune201716S
## 1093 1089 1088
## 45RRNMarch201716S 45RRNMarch201816S 45RRNSeptember201716S
## 1090 1092 1091
## 46LXNAugust201616S 46LXNDecember201616S 46LXNDecember201716S
## 1089 1090 1092
## 46LXNJune201716S 46LXNMarch201716S 46LXNMarch201816S
## 1087 1093 1090
## 46RXXAugust201616S 46RXXDecember201616S 46RXXJune201716S
## 1092 1091 1091
## 46RXXMarch201716S 46RXXMarch201816S 46RXXSeptember201716S
## 1089 1092 1090
## 47LRNAugust201616S 47LRNDecember201616S 47LRNJune201716S
## 1090 1088 1090
## 47LRNMarch201716S 47LRNMarch201816S 47LRNSeptember201716S
## 1091 1090 1090
## 47RRXAugust201616S 47RRXDecember201616S 47RRXJune201716S
## 1090 1091 1093
## 47RRXMarch201716S 47RRXMarch201816S 47RRXSeptember201716S
## 1090 1086 1090
## 48LXXAugust201616S 48LXXDecember201616S 48LXXDecember201716S
## 1089 1088 1087
## 48LXXJune201716S 48LXXMarch201716S 48LXXMarch201816S
## 1088 1089 1088
## 48RXNAugust201616S 48RXNDecember201616S 48RXNJune201716S
## 1087 1086 1090
## 48RXNMarch201716S 48RXNMarch201816S 48RXNSeptember201716S
## 1092 1091 1092
## 4LXXDecember201616S 4LXXDecember201716S 4LXXMarch201716S
## 1090 1090 1089
## 4LXXMarch201816S 4RXNAugust201616S 4RXNDecember201616S
## 1089 1088 1089
## 4RXNDecember201716S 4RXNJune201716S 4RXNMarch201716S
## 1089 1090 1087
## 4RXNMarch201816S 4RXNSeptember201716S 5LRNDecember201616S
## 1090 1089 1091
## 5LRNDecember201716S 5LRNMarch201716S 5LRNMarch201816S
## 1090 1089 1091
## 5LRNSeptember201716S 5RRXAugust201616S 5RRXDecember201616S
## 1090 1089 1091
## 5RRXJune201716S 5RRXMarch201716S 5RRXMarch201816S
## 1089 1089 1091
## 5RRXSeptember201716S 7LRXDecember201616S 7LRXDecember201716S
## 1091 1093 1090
## 7LRXMarch201716S 7LRXMarch201816S 7LRXSeptember201716S
## 1092 1090 1090
## 7RRNDecember201616S 7RRNDecember201716S 7RRNJune201716S
## 1091 1089 1091
## 7RRNMarch201716S 7RRNMarch201816S 7RRNSeptember201716S
## 1092 1091 1090
## 8LXXAugust201616S 8LXXDecember201616S 8LXXDecember201716S
## 1091 1089 1091
## 8LXXMarch201716S 8LXXMarch201816S 8LXXSeptember201716S
## 1090 1088 1090
## 8RXNAugust201616S 8RXNDecember201616S
## 1088 1090
# Check the number of unique OTUs at different levels
phylum.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Phylum[!is.na(otu.raref.taxonomy.16s.ss$Phylum)]))
length(phylum.16s) #15 phyla
## [1] 15
class.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Class[!is.na(otu.raref.taxonomy.16s.ss$Class)]))
length(class.16s) #36 unique bacterial classes (NA not counted)
## [1] 34
order.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Order[!is.na(otu.raref.taxonomy.16s.ss$Order)]))
length(order.16s) #55 unique bacterial orders (na not counted)
## [1] 53
family.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Family[!is.na(otu.raref.taxonomy.16s.ss$Family)]))
length(family.16s) #97 unique families (na not counted)
## [1] 96
genus.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Genus[!is.na(otu.raref.taxonomy.16s.ss$Genus)]))
length(genus.16s) #142 unique bacterial genus' (na not counted)
## [1] 141
#aggregate taxonomical IDs on a certain level by sample (here at Order level)
genus.16s.1 <- aggregate(. ~Genus, data = otu.raref.taxonomy.16s.ss[c(1:(length(names(otu.raref.taxonomy.16s.ss)) - 7), match("Genus", names(otu.raref.taxonomy.16s.ss)))], FUN = sum, na.omit=TRUE)
#Order column is turned into rownames and the empty cell is called "unidentified" (o__ gsub into "")
rownames(genus.16s.1) <- genus.16s.1$Genus
genus.16s.1$Genus <- NULL
#turn table into dataframe and transpose
genus.16s.1 <- as.data.frame(t(genus.16s.1))
write.xlsx(genus.16s.1, file = "genera.16s.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names=TRUE)
# Merge a metadata column here (I just made my own metatdata column)
matching.metadata.gen.16s <- metadata.16s.1[c(match(sort(row.names(genus.16s.1)), sort(row.names(metadata.16s.1)))),]
merged.meta.gen.16s <- merge(genus.16s.1, matching.metadata.gen.16s, by.x = "row.names", by.y = "row.names")
colnames(merged.meta.gen.16s)[colnames(merged.meta.gen.16s)=="Row.names"] <- "SampleID"
merged.meta.gen.16s <- as.data.frame(merged.meta.gen.16s)
rownames(merged.meta.gen.16s) <- merged.meta.gen.16s$SampleID
merged.meta.gen.16s$SampleID <- NULL
#aggregate only the X number of the level of taxonomy chosen (from the sequence summaries above) with the metadata selected in this case plot Treatment
genera.agg.16s <- aggregate(merged.meta.gen.16s[, 1:142], list(merged.meta.gen.16s$eco.prec), mean) #genus
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
genera.agg.16s <- as.data.frame(t(genera.agg.16s)) # transpose table to have variables as columns
names(genera.agg.16s) <- as.character(unlist(genera.agg.16s["Group.1", ])) # transposing causes the treatments named as "Group.1", remove this
genera.agg.16s <- genera.agg.16s[-1,] #final removal of "Group.1"
for(i in c(1:ncol(genera.agg.16s))) {
genera.agg.16s[,i] <- as.numeric(as.character(genera.agg.16s[,i]))
}
write.xlsx(genera.agg.16s, file = "plot.genera.16s.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names=TRUE)
#calculate the number of sequences per level (change the divisor according to above calculation for unique OTUs per level)
totalseqs.16s <- colSums(genera.agg.16s)
avgnumseqs.16s <- sum(totalseqs.16s)/142 #142 is the total number of unique genera with NA removed
# CHECK FOR THE MOST ABUNDANT of classifaction level selected HERE!!!!
genera.agg.16s$Genus <- row.names(genera.agg.16s)
genera.agg.16s$Genus <- factor(genera.agg.16s$Genus)
sumsum.16s <- colSums(genera.agg.16s[1:(length(names(genera.agg.16s))-1)])
samp.16s <- names(genera.agg.16s)[1:(length(names(genera.agg.16s))-1)]
genera.agg.long.16s <- melt(genera.agg.16s, id.vars = "Genus", measure.vars = c("CSS_Ambient","CSS_Reduced", "Grass_Ambient","Grass_Reduced"))
## Warning in melt(genera.agg.16s, id.vars = "Genus", measure.vars =
## c("CSS_Ambient", : The melt generic in data.table has been passed a data.frame
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(genera.agg.
## 16s). In the next version, this warning will become an error.
genera.agg.long.16s$Abundance <- ifelse(genera.agg.long.16s$variable == samp.16s[1], genera.agg.long.16s$value / sumsum.16s[1],
ifelse(genera.agg.long.16s$variable == samp.16s[2], genera.agg.long.16s$value / sumsum.16s[2],
ifelse(genera.agg.long.16s$variable == samp.16s[3], genera.agg.long.16s$value / sumsum.16s[3],
ifelse(genera.agg.long.16s$variable == samp.16s[4], genera.agg.long.16s$value / sumsum.16s[4],"NA" ))))
#call everything below 0.01 (1%) "other_genera"
genera.agg.long.16s$Genus[genera.agg.long.16s$Abundance < 0.01] <- genera.agg.long.16s$Genus["other_genera"] #turn all genera that are below 1% first into "NA" then into "other_genera"
genera.agg.long.16s$Genus <- as.character(genera.agg.long.16s$Genus)
genera.agg.long.16s$Genus[is.na(genera.agg.long.16s$Genus)] <- "other_genera"
genera.agg.long.16s$Abundance <- as.numeric(genera.agg.long.16s$Abundance)
# Checking to make sure it adds to 1 precisely
for (i in 1:length(samp.16s)) {
print(paste0("For precipitation treatment ", samp.16s[i], " the relative abundance adds to: ", sum(genera.agg.long.16s[genera.agg.long.16s$variable == samp.16s[i], ]$Abundance), "."))
}
## [1] "For precipitation treatment CSS_Ambient the relative abundance adds to: NA."
## [1] "For precipitation treatment CSS_Reduced the relative abundance adds to: NA."
## [1] "For precipitation treatment Grass_Ambient the relative abundance adds to: NA."
## [1] "For precipitation treatment Grass_Reduced the relative abundance adds to: NA."
#create stacked barplot
genera.agg.long.16s$Genus <- factor(genera.agg.long.16s$Genus)
genera.agg.long.16s$Genus = with(genera.agg.long.16s, reorder(Genus, -Abundance))
#genera.agg.long.16s$Genus <- relevel(genera.agg.long.16s$Genus, 'other_genera') #if want the "other_genera" at the top of the plot.
c36 <- c("#c8c9c9",
"#bda69d",
"#284d50",
"#c1a4a6",
"#344b4a",
"#b2a6bd",
"#58423a",
"#8eb2a9",
"#594152",
"#83adad",
"#56414b",
"#b2aa97",
"#484555",
"#a9ac9f",
"#3a4951",
"#b4a8ab",
"#414849",
"#a1adac",
"#544341",
"#749ca7",
"#5f5740",
"#a98d9e",
"#4b5d49",
"#b4908d",
"#3a586a",
"#896d5f",
"#4e5269",
"#7d8168",
"#665b72",
"#718c7a",
"#6d505d",
"#517368",
"#86656c",
"#46474a",
"#7c819a",
"#967471")
barplot.genus.16s <- ggplot(genera.agg.long.16s, aes(x = fct_relevel(variable, "Grass_Ambient","Grass_Reduced","CSS_Ambient","CSS_Reduced"), y = Abundance, fill = Genus)) +
theme_test() +
geom_bar(stat = "identity", width = 0.7) +
theme(axis.line.x=element_line(colour = "black"), axis.line.y=element_line(colour = "black")) +
scale_fill_manual(values = c36) +
scale_y_continuous(expand = c(0, 0)) +
labs(fill = "Genus") +
theme(legend.position = "none") + #Change to "right" to display legend at the right of the plot
guides(size = "none") +
labs(x = "Ecosystem by Water Input", y = "Proportion of Total Bacterial Community")
plot(barplot.genus.16s)
## Warning: Removed 568 rows containing missing values (position_stack).

###Fungi###
taxonomyITS<-read_qza("taxonomyITS.qza", "./tmp", rm=T)
taxonomyITS$uuid
## [1] "e4992fa1-ce2d-4474-b859-f5103b8c08e2"
taxtableITS<-taxonomyITS$data %>% as.tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species")) #convert the table into a tabular split version
## Warning: Expected 7 pieces. Missing pieces filled with `NA` in 2197 rows [1, 3,
## 4, 5, 6, 7, 8, 10, 11, 13, 16, 17, 20, 21, 22, 23, 24, 25, 26, 27, ...].
taxtableITS
## # A tibble: 4,601 x 9
## Feature.ID Kingdom Phylum Class Order Family Genus Species Confidence
## <fct> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 7e80c7e5f58… k__Fungi p__Asc… <NA> <NA> <NA> <NA> <NA> 0.812
## 2 1094975b0bd… k__Fungi p__Asc… c__Le… o__H… f__Hel… g__Ar… s__Arti… 0.952
## 3 3561d8e1d04… k__Fungi p__Asc… c__So… o__G… f__Glo… g__Co… <NA> 0.982
## 4 9995cbaaec4… k__Fungi p__Asc… <NA> <NA> <NA> <NA> <NA> 0.820
## 5 592482291be… k__Fungi p__Asc… <NA> <NA> <NA> <NA> <NA> 0.836
## 6 7d8e2be52cd… k__Fungi p__Asc… c__Do… o__C… <NA> <NA> <NA> 0.982
## 7 0ca5bba7ce1… k__Fungi p__Asc… c__Do… <NA> <NA> <NA> <NA> 0.836
## 8 53f1a5397d6… k__Plan… p__Chl… c__Tr… o__T… f__Tre… g__Tr… <NA> 0.999
## 9 210acd46015… k__Fungi p__Asc… c__Do… o__P… f__Pha… g__Ch… s__unid… 0.995
## 10 fb069049f2a… k__Fungi <NA> <NA> <NA> <NA> <NA> <NA> 0.984
## # … with 4,591 more rows
taxtableITS <- as.data.frame(taxtableITS)
#remove p__ as a precursor to phylum name
taxtableITS$Phylum <- gsub("p__", "", taxtableITS$Phylum)
taxtableITS$Class <- gsub("c__", "", taxtableITS$Class)
taxtableITS$Order <- gsub("o__", "", taxtableITS$Order)
taxtableITS$Family <- gsub("f__", "", taxtableITS$Family)
taxtableITS$Genus <- gsub("g__", "", taxtableITS$Genus)
taxtableITS$Species <- gsub("s__", "", taxtableITS$Species)
#this will turn all empty cell into NA
taxtableITS <- taxtableITS %>%
mutate_at(vars(colnames(.)),
.funs = funs(ifelse(.=="", NA, as.character(.))))
#taxtable was imported using qiime2R, with already separated levels
row.names(taxtableITS) <- taxtableITS$Feature.ID #turn feature IDs into row.names
taxtableITS$Feature.ID <- NULL #remove column "Feature ID"
#merge OTU table with taxonomy table
otu.raref.taxonomy.ITS <- as.data.frame(merge(taxtableITS, rared.otu.ITS.t, by.x = "row.names", by.y = "row.names"))
row.names(otu.raref.taxonomy.ITS) <- otu.raref.taxonomy.ITS$Row.names #move column "Row.names" (first column) to row.names this was created during the merge
otu.raref.taxonomy.ITS$Row.names <- NULL #remove column "Row.names"
otu.raref.taxonomy.ITS$Confidence <- NULL # remove column "Confidence"
otu.raref.taxonomy.ITS.ss <- subset(otu.raref.taxonomy.ITS, select=c(8:195,1:7)) # move Taxa information to the end of all columns
write.table(otu.raref.taxonomy.ITS.ss, "otu_raref_taxonomy.ITS.csv", quote=FALSE, sep=',') #use this table to calculate top number of ESVs and their sequence contributions.
colSums(otu.raref.taxonomy.ITS.ss[, c(1:188)]) # check if samples are RAREFIED!!
## 14LRNDecember2016ITS 14LRNDecember2017ITS 14LRNJune2017ITS
## 1066 1063 1066
## 14LRNMarch2017ITS 14LRNMarch2018ITS 14LRNSeptember2017ITS
## 1065 1064 1068
## 14RRXAugust2016ITS 14RRXDecember2016ITS 14RRXDecember2017ITS
## 1067 1060 1068
## 14RRXJune2017ITS 14RRXMarch2017ITS 14RRXMarch2018ITS
## 1063 1066 1067
## 14RRXSeptember2017ITS 18LXNAugust2016ITS 18LXNDecember2016ITS
## 1063 1065 1062
## 18LXNDecember2017ITS 18LXNJune2017ITS 18LXNMarch2017ITS
## 1065 1061 1068
## 18LXNMarch2018ITS 18LXNSeptember2017ITS 18RXXAugust2016ITS
## 1064 1063 1063
## 18RXXDecember2016ITS 18RXXDecember2017ITS 18RXXJune2017ITS
## 1064 1068 1064
## 18RXXMarch2017ITS 18RXXMarch2018ITS 18RXXSeptember2017ITS
## 1067 1058 1064
## 20LRXAugust2016ITS 20LRXDecember2016ITS 20LRXDecember2017ITS
## 1065 1063 1062
## 20LRXJune2017ITS 20LRXMarch2017ITS 20LRXMarch2018ITS
## 1063 1064 1063
## 20LRXSeptember2017ITS 20RRNDecember2016ITS 20RRNDecember2017ITS
## 1065 1068 1066
## 20RRNMarch2017ITS 20RRNMarch2018ITS 20RRNSeptember2017ITS
## 1067 1067 1066
## 22LXXAugust2016ITS 22LXXDecember2016ITS 22LXXDecember2017ITS
## 1068 1068 1064
## 22LXXJune2017ITS 22LXXMarch2017ITS 22LXXMarch2018ITS
## 1063 1066 1063
## 22LXXSeptember2017ITS 22RXNAugust2016ITS 22RXNDecember2016ITS
## 1064 1061 1061
## 22RXNDecember2017ITS 22RXNJune2017ITS 22RXNMarch2017ITS
## 1062 1066 1063
## 22RXNMarch2018ITS 22RXNSeptember2017ITS 25LRXAugust2016ITS
## 1063 1064 1066
## 25LRXDecember2016ITS 25LRXDecember2017ITS 25LRXMarch2017ITS
## 1061 1067 1064
## 25LRXMarch2018ITS 25LRXSeptember2017ITS 25RRNAugust2016ITS
## 1060 1061 1066
## 25RRNDecember2016ITS 25RRNDecember2017ITS 25RRNJune2017ITS
## 1061 1063 1065
## 25RRNMarch2017ITS 25RRNMarch2018ITS 25RRNSeptember2017ITS
## 1059 1066 1065
## 27LXNAugust2016ITS 27LXNDecember2016ITS 27LXNDecember2017ITS
## 1067 1064 1065
## 27LXNMarch2017ITS 27LXNMarch2018ITS 27LXNSeptember2017ITS
## 1066 1063 1062
## 27RXXDecember2016ITS 27RXXJune2017ITS 27RXXMarch2017ITS
## 1062 1062 1064
## 27RXXMarch2018ITS 32LXNAugust2016ITS 32LXNDecember2016ITS
## 1065 1064 1065
## 32LXNMarch2017ITS 32LXNMarch2018ITS 32LXNSeptember2017ITS
## 1061 1062 1061
## 32RXXAugust2016ITS 32RXXDecember2016ITS 32RXXDecember2017ITS
## 1067 1063 1064
## 32RXXJune2017ITS 32RXXMarch2017ITS 32RXXMarch2018ITS
## 1068 1059 1066
## 32RXXSeptember2017ITS 35LRXAugust2016ITS 35LRXDecember2016ITS
## 1058 1064 1061
## 35LRXDecember2017ITS 35LRXJune2017ITS 35LRXMarch2017ITS
## 1067 1067 1066
## 35LRXMarch2018ITS 35LRXSeptember2017ITS 35RRNAugust2016ITS
## 1065 1064 1059
## 35RRNDecember2016ITS 35RRNDecember2017ITS 35RRNJune2017ITS
## 1060 1066 1059
## 35RRNMarch2017ITS 35RRNMarch2018ITS 35RRNSeptember2017ITS
## 1061 1065 1059
## 45LRXDecember2016ITS 45LRXDecember2017ITS 45LRXJune2017ITS
## 1064 1066 1062
## 45LRXMarch2017ITS 45LRXMarch2018ITS 45LRXSeptember2017ITS
## 1062 1065 1062
## 45RRNAugust2016ITS 45RRNDecember2016ITS 45RRNDecember2017ITS
## 1067 1066 1064
## 45RRNJune2017ITS 45RRNMarch2017ITS 45RRNMarch2018ITS
## 1061 1067 1062
## 45RRNSeptember2017ITS 46LXNAugust2016ITS 46LXNDecember2016ITS
## 1067 1059 1063
## 46LXNDecember2017ITS 46LXNJune2017ITS 46LXNMarch2017ITS
## 1065 1069 1062
## 46LXNMarch2018ITS 46LXNSeptember2017ITS 46RXXAugust2016ITS
## 1067 1063 1065
## 46RXXDecember2016ITS 46RXXDecember2017ITS 46RXXJune2017ITS
## 1067 1065 1061
## 46RXXMarch2018ITS 46RXXSeptember2017ITS 47LRNAugust2016ITS
## 1063 1067 1066
## 47LRNDecember2016ITS 47LRNDecember2017ITS 47LRNJune2017ITS
## 1064 1064 1063
## 47LRNMarch2017ITS 47LRNMarch2018ITS 47LRNSeptember2017ITS
## 1064 1068 1065
## 47RRXAugust2016ITS 47RRXDecember2016ITS 47RRXJune2017ITS
## 1064 1063 1061
## 47RRXMarch2017ITS 47RRXMarch2018ITS 47RRXSeptember2017ITS
## 1063 1066 1067
## 48LXXAugust2016ITS 48LXXDecember2016ITS 48LXXDecember2017ITS
## 1068 1064 1062
## 48LXXJune2017ITS 48LXXMarch2017ITS 48LXXMarch2018ITS
## 1067 1066 1072
## 48RXNAugust2016ITS 48RXNDecember2016ITS 48RXNDecember2017ITS
## 1066 1065 1068
## 48RXNJune2017ITS 48RXNMarch2017ITS 48RXNMarch2018ITS
## 1065 1058 1071
## 48RXNSeptember2017ITS 4LXXAugust2016ITS 4LXXDecember2016ITS
## 1067 1064 1066
## 4LXXDecember2017ITS 4LXXJune2017ITS 4LXXMarch2017ITS
## 1064 1061 1066
## 4LXXSeptember2017ITS 4RXNAugust2016ITS 4RXNDecember2017ITS
## 1059 1067 1064
## 4RXNJune2017ITS 4RXNMarch2017ITS 4RXNMarch2018ITS
## 1063 1063 1065
## 4RXNSeptember2017ITS 5LRNDecember2016ITS 5LRNJune2017ITS
## 1063 1067 1062
## 5LRNMarch2017ITS 5LRNMarch2018ITS 5LRNSeptember2017ITS
## 1059 1064 1061
## 5RRXAugust2016ITS 5RRXDecember2016ITS 5RRXDecember2017ITS
## 1057 1065 1064
## 5RRXJune2017ITS 5RRXMarch2017ITS 5RRXMarch2018ITS
## 1071 1064 1066
## 5RRXSeptember2017ITS 7LRXAugust2016ITS 7LRXDecember2016ITS
## 1066 1061 1065
## 7LRXDecember2017ITS 7LRXJune2017ITS 7LRXMarch2017ITS
## 1066 1067 1064
## 7LRXMarch2018ITS 7LRXSeptember2017ITS 7RRNDecember2016ITS
## 1067 1063 1066
## 7RRNDecember2017ITS 7RRNJune2017ITS
## 1060 1067
# Check the number of unique OTUs at different levels
phylum.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Phylum[!is.na(otu.raref.taxonomy.ITS.ss$Phylum)]))
length(phylum.ITS) #4 phyla
## [1] 4
class.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Class[!is.na(otu.raref.taxonomy.ITS.ss$Class)]))
length(class.ITS) #21 unique bacterial classes (NA not counted)
## [1] 21
order.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Order[!is.na(otu.raref.taxonomy.ITS.ss$Order)]))
length(order.ITS) #67 unique bacterial orders (na not counted)
## [1] 67
family.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Family[!is.na(otu.raref.taxonomy.ITS.ss$Family)]))
length(family.ITS) #155 unique families (na not counted)
## [1] 155
genus.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Genus[!is.na(otu.raref.taxonomy.ITS.ss$Genus)]))
length(genus.ITS) #275 unique bacterial genus' (na not counted)
## [1] 274
#aggregate taxonomical IDs on a certain level by sample (here at Order level)
genus.ITS.1 <- aggregate(. ~Genus, data = otu.raref.taxonomy.ITS.ss[c(1:(length(names(otu.raref.taxonomy.ITS.ss)) - 7), match("Genus", names(otu.raref.taxonomy.ITS.ss)))], FUN = sum, na.omit=TRUE)
#Order column is turned into rownames and the empty cell is called "unidentified" (o__ gsub into "")
rownames(genus.ITS.1) <- genus.ITS.1$Genus
genus.ITS.1$Genus <- NULL
#turn table into dataframe and transpose
genus.ITS.1 <- as.data.frame(t(genus.ITS.1))
write.xlsx(genus.ITS.1, file = "genera.ITS.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names=TRUE)
# Merge a metadata column here (I just made my own metatdata column)
matching.metadata.gen.ITS <- metadata.ITS.1[c(match(sort(row.names(genus.ITS.1)), sort(row.names(metadata.ITS.1)))),]
merged.meta.gen.ITS <- merge(genus.ITS.1, matching.metadata.gen.ITS, by.x = "row.names", by.y = "row.names")
colnames(merged.meta.gen.ITS)[colnames(merged.meta.gen.ITS)=="Row.names"] <- "SampleID"
merged.meta.gen.ITS <- as.data.frame(merged.meta.gen.ITS)
rownames(merged.meta.gen.ITS) <- merged.meta.gen.ITS$SampleID
merged.meta.gen.ITS$SampleID <- NULL
#aggregate only the X number of the level of taxonomy chosen (from the sequence summaries above) with the metadata selected in this case plot Treatment
genera.agg.ITS <- aggregate(merged.meta.gen.ITS[, 1:275], list(merged.meta.gen.ITS$eco.prec), mean) #genus
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
genera.agg.ITS <- as.data.frame(t(genera.agg.ITS)) # transpose table to have variables as columns
names(genera.agg.ITS) <- as.character(unlist(genera.agg.ITS["Group.1", ])) # transposing causes the treatments named as "Group.1", remove this
genera.agg.ITS <- genera.agg.ITS[-1,] #final removal of "Group.1"
for(i in c(1:ncol(genera.agg.ITS))) {
genera.agg.ITS[,i] <- as.numeric(as.character(genera.agg.ITS[,i]))
}
write.xlsx(genera.agg.ITS, file = "plot.genera.ITS.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names=TRUE)
#calculate the number of sequences per level (change the divisor according to above calculation for unique OTUs per level)
totalseqs.ITS <- colSums(genera.agg.ITS)
avgnumseqs.ITS <- sum(totalseqs.ITS)/275 #142 is the total number of unique genera with NA removed
# CHECK FOR THE MOST ABUNDANT of classifaction level selected HERE!!!!
genera.agg.ITS$Genus <- row.names(genera.agg.ITS)
genera.agg.ITS$Genus <- factor(genera.agg.ITS$Genus)
sumsum.ITS <- colSums(genera.agg.ITS[1:(length(names(genera.agg.ITS))-1)])
samp.ITS <- names(genera.agg.ITS)[1:(length(names(genera.agg.ITS))-1)]
genera.agg.long.ITS <- melt(genera.agg.ITS, id.vars = "Genus", measure.vars = c("CSS_Ambient","CSS_Reduced", "Grass_Ambient","Grass_Reduced"))
## Warning in melt(genera.agg.ITS, id.vars = "Genus", measure.vars =
## c("CSS_Ambient", : The melt generic in data.table has been passed a
## data.frame and will attempt to redirect to the relevant reshape2 method;
## please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(genera.agg.ITS). In the next version, this warning will become an
## error.
genera.agg.long.ITS$Abundance <- ifelse(genera.agg.long.ITS$variable == samp.ITS[1], genera.agg.long.ITS$value / sumsum.ITS[1],
ifelse(genera.agg.long.ITS$variable == samp.ITS[2], genera.agg.long.ITS$value / sumsum.ITS[2],
ifelse(genera.agg.long.ITS$variable == samp.ITS[3], genera.agg.long.ITS$value / sumsum.ITS[3],
ifelse(genera.agg.long.ITS$variable == samp.ITS[4], genera.agg.long.ITS$value / sumsum.ITS[4],"NA" ))))
#call everything below 0.01 (1%) "other_genera"
genera.agg.long.ITS$Genus[genera.agg.long.ITS$Abundance < 0.01] <- genera.agg.long.ITS$Genus["other_genera"] #turn all genera that are below 1% first into "NA" then into "other_genera"
genera.agg.long.ITS$Genus[ genera.agg.long.ITS$Genus == "unidentified" ] <- NA
genera.agg.long.ITS$Genus <- as.character(genera.agg.long.ITS$Genus)
genera.agg.long.ITS$Genus[is.na(genera.agg.long.ITS$Genus)] <- "other_genera"
genera.agg.long.ITS$Abundance <- as.numeric(genera.agg.long.ITS$Abundance)
# Checking to make sure it adds to 1 precisely
for (i in 1:length(samp.ITS)) {
print(paste0("For precipitation treatment ", samp.ITS[i], " the relative abundance adds to: ", sum(genera.agg.long.ITS[genera.agg.long.ITS$variable == samp.ITS[i], ]$Abundance), "."))
}
## [1] "For precipitation treatment CSS_Ambient the relative abundance adds to: NA."
## [1] "For precipitation treatment CSS_Reduced the relative abundance adds to: NA."
## [1] "For precipitation treatment Grass_Ambient the relative abundance adds to: NA."
## [1] "For precipitation treatment Grass_Reduced the relative abundance adds to: NA."
#create stacked barplot
genera.agg.long.ITS$Genus <- factor(genera.agg.long.ITS$Genus)
genera.agg.long.ITS$Genus = with(genera.agg.long.ITS, reorder(Genus, -Abundance))
#genera.agg.long.ITS$Genus <- relevel(genera.agg.long.ITS$Genus, "other_genera") #if want the "other_genera" at the top of the plot.
c18 <- c("#c8c9c9",
"#264d4c",
"#bfa4ab",
"#88abbc",
"#645b44",
"#ac9fb8",
"#736f56",
"#9692ac",
"#7c6054",
"#5f8791",
"#b7988b",
"#604b5e",
"#738b77",
"#93737e",
"#5a7c71",
"#8e8c72",
"#587789",
"#676666")
barplot.genus.ITS <- ggplot(genera.agg.long.ITS, aes(x = fct_relevel(variable, "Grass_Ambient","Grass_Reduced","CSS_Ambient","CSS_Reduced"), y = Abundance, fill = Genus)) +
theme_test() +
geom_bar(stat = "identity", width = 0.7) +
theme(axis.line.x=element_line(colour = "black"), axis.line.y=element_line(colour = "black")) +
scale_fill_manual(values = c18) +
scale_y_continuous(expand = c(0, 0)) +
labs(fill = "Genus") +
theme(legend.position = "none") + #Change to "right" to display legend at the right of the plot
guides(size = "none") +
labs(x = "Ecosystem by Water Input", y = "Proportion of Total Fungal Community")
plot(barplot.genus.ITS)
## Warning: Removed 1100 rows containing missing values (position_stack).

###Plant###
plant.func.2015 <- read_excel("plant_barplot.2015.eco.pre.xlsx")
plant.func.2015 <- as.data.frame(plant.func.2015, col.names = TRUE, row.names = TRUE)
plant.func.2015[1:3 , 1:3]
## functional_group community percentage
## 1 g_Bare css_ambient 2.988300
## 2 g_Bare css_reduced 25.266378
## 3 g_Bare grass_ambient 2.143758
c07 <- c("#344c60",
"#baa6ab",
"#3f4944",
"#7b9a8b",
"#85778e",
"#5c5f47",
"#82655a")
barplot.plant.2015 <- ggplot(plant.func.2015, aes(x = fct_relevel(community,"grass_ambient","grass_reduced","css_ambient","css_reduced"), y = percentage, fill = functional_group)) +
geom_col(aes(fill = factor(functional_group, levels = c("NativeGrassCover", "NonNativeGrassCover","NativeShrubCover","NativeForbCover", "NonNativeForbCover", "Litter",
"Bare")))) +
theme_test() +
geom_bar(stat = "identity", width = 0.7) +
theme(axis.line.x=element_line(colour = "black"), axis.line.y=element_line(colour = "black")) +
scale_y_continuous(expand = c(0, 0)) +
scale_fill_manual(values = c07) +
labs(fill = "factor") +
theme(legend.position = "none") + #Change to "right" to display legend at the right of the plot
guides(size = "right") +
labs(x = "Plant Functional Groups by Ecosystem", y = "Proportion Fractional Cover")
plot(barplot.plant.2015)

###Plot Bacterial, Fungal, and Plant in one figure###
combined.stacked.taxa.genus <- ggarrange(barplot.genus.16s, barplot.genus.ITS, barplot.plant.2015, barplot.plant.2015, labels = c("A.", "B.","C.","D."), ncol = 2, nrow = 2, widths = c(5,5,5,5))
## Warning: Removed 568 rows containing missing values (position_stack).
## Warning: Removed 1100 rows containing missing values (position_stack).
plot(combined.stacked.taxa.genus)

###Richness and Evenness###
###Bacteria###
attach(metadata.16s.1)
Nitrogen<-factor(Nitrogen)
Precipitation <- as.factor(Precipitation)
Block <- as.factor(Block)
Plot <- as.factor(Plot)
Collection_date <- as.factor(Collection_date)
eco.prec <- as.factor(eco.prec)
eco.nit <- as.factor(eco.nit)
bact.aov.richness <- lmer(richness ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen + ( 1 | Block:Ecosystem), REML = FALSE, data = metadata.16s.1)
## boundary (singular) fit: see ?isSingular
anova(bact.aov.richness, type = "III")
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Block 2279.8 2279.8 1 192 7.7321 0.0059650 **
## Ecosystem 3545.2 3545.2 1 192 12.0241 0.0006481 ***
## Precipitation 74.4 74.4 1 192 0.2523 0.6160429
## Nitrogen 78.8 78.8 1 192 0.2674 0.6056878
## Collection_date 12177.5 2029.6 6 192 6.8837 1.234e-06 ***
## Ecosystem:Precipitation 526.7 526.7 1 192 1.7863 0.1829637
## Ecosystem:Nitrogen 91.9 91.9 1 192 0.3117 0.5772818
## Ecosystem:Collection_date 2552.8 425.5 6 192 1.4430 0.2001184
## Precipitation:Collection_date 3396.3 566.1 6 192 1.9199 0.0794684 .
## Nitrogen:Collection_date 1137.2 189.5 6 192 0.6428 0.6958267
## Precipitation:Nitrogen 1373.4 1373.4 1 192 4.6580 0.0321501 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bact.eco.prec.rich <- ggplot(data = bact.aov.richness, aes(x = factor(eco.prec, levels = c("Grass_Ambient","Grass_Reduced","CSS_Ambient", "CSS_Reduced")), y = metadata.16s.1$richness, fill = metadata.16s.1$Precipitation)) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center',
position=position_dodge(0.75), dotsize = 0.65) +
scale_fill_manual(values=c('#408080', '#FF8040')) +
labs(title = 'Bacterial Richness by Ecosystem', x = 'Ecosystem', y = 'Richness', fill = 'Precipitation') +
theme_classic() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_blank())
plot(bact.eco.prec.rich)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

bact.eco.prec.tuk.rich <- TukeyHSD(aov(formula = metadata.16s.1$richness ~ as.character(metadata.16s.1$eco.prec)))
bact.eco.prec.tuk.rich
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = metadata.16s.1$richness ~ as.character(metadata.16s.1$eco.prec))
##
## $`as.character(metadata.16s.1$eco.prec)`
## diff lwr upr p adj
## CSS_Reduced-CSS_Ambient -4.9400000 -15.657512 5.777512 0.6308928
## Grass_Ambient-CSS_Ambient -7.7208333 -18.549410 3.107744 0.2542266
## Grass_Reduced-CSS_Ambient -4.8363636 -15.913221 6.240494 0.6703734
## Grass_Ambient-CSS_Reduced -2.7808333 -13.609410 8.047744 0.9098615
## Grass_Reduced-CSS_Reduced 0.1036364 -10.973221 11.180494 0.9999949
## Grass_Reduced-Grass_Ambient 2.8844697 -8.299885 14.068824 0.9088281
bact.eco.nit.rich <- ggplot(data = bact.aov.richness, aes(x = factor(eco.nit, levels = c("Grass_Ambient","Grass_Added","CSS_Ambient", "CSS_Added")), y = metadata.16s.1$richness, fill = metadata.16s.1$Nitrogen)) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center',
position=position_dodge(0.75), dotsize = 0.65) +
scale_fill_manual(values=c('#b7a28c','#408080')) +
labs(title = 'Bacterial Richness by Ecosystem', x = 'Ecosystem', y = 'Richness', fill = 'Nitrogen') +
theme_classic() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_blank())
plot(bact.eco.nit.rich)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

bact.aov.shannon <- lmer(shannon ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen + ( 1 | Block:Ecosystem), REML = FALSE, data = metadata.16s.1)
anova(bact.aov.shannon, type = "III")
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Block 0.8067 0.80671 1 7.682 3.5760 0.096799
## Ecosystem 1.9936 1.99360 1 7.609 8.8373 0.018832
## Precipitation 0.2835 0.28350 1 184.221 1.2567 0.263732
## Nitrogen 0.1202 0.12023 1 186.131 0.5330 0.466284
## Collection_date 9.5941 1.59902 6 184.997 7.0882 8.248e-07
## Ecosystem:Precipitation 0.0268 0.02678 1 184.032 0.1187 0.730827
## Ecosystem:Nitrogen 0.1979 0.19785 1 185.400 0.8771 0.350228
## Ecosystem:Collection_date 2.2662 0.37771 6 185.186 1.6743 0.129466
## Precipitation:Collection_date 2.9398 0.48997 6 184.617 2.1719 0.047575
## Nitrogen:Collection_date 1.2565 0.20941 6 185.166 0.9283 0.475847
## Precipitation:Nitrogen 1.8126 1.81260 1 184.929 8.0349 0.005099
##
## Block .
## Ecosystem *
## Precipitation
## Nitrogen
## Collection_date ***
## Ecosystem:Precipitation
## Ecosystem:Nitrogen
## Ecosystem:Collection_date
## Precipitation:Collection_date *
## Nitrogen:Collection_date
## Precipitation:Nitrogen **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bact.eco.prec.shannon <- ggplot(data = bact.aov.shannon, aes(x = factor(eco.prec, levels = c("Grass_Ambient","Grass_Reduced","CSS_Ambient", "CSS_Reduced")), y = metadata.16s.1$shannon, fill = metadata.16s.1$Precipitation)) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center',
position=position_dodge(0.75), dotsize = 0.65) +
scale_fill_manual(values=c('#408080', '#FF8040')) +
labs(title = 'Bacterial Evenness by Ecosystem', x = 'Ecosystem', y = 'Evenness', fill = 'Precipitation') +
theme_classic() +
#theme(legend.position = "none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())
bact.eco.prec.tuk.shannon <- TukeyHSD(aov(formula = metadata.16s.1$shannon ~ as.character(metadata.16s.1$eco.prec)))
bact.eco.prec.tuk.shannon
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = metadata.16s.1$shannon ~ as.character(metadata.16s.1$eco.prec))
##
## $`as.character(metadata.16s.1$eco.prec)`
## diff lwr upr p adj
## CSS_Reduced-CSS_Ambient -0.06079776 -0.3615575 0.23996197 0.9532127
## Grass_Ambient-CSS_Ambient -0.16421605 -0.4680925 0.13966044 0.5004130
## Grass_Reduced-CSS_Ambient -0.23491663 -0.5457605 0.07592721 0.2073915
## Grass_Ambient-CSS_Reduced -0.10341829 -0.4072948 0.20045820 0.8141064
## Grass_Reduced-CSS_Reduced -0.17411887 -0.4849627 0.13672497 0.4687036
## Grass_Reduced-Grass_Ambient -0.07070058 -0.3845611 0.24315991 0.9368327
bact.eco.nit.shannon <- ggplot(data = bact.aov.shannon, aes(x = factor(eco.nit, levels = c("Grass_Ambient","Grass_Added","CSS_Ambient", "CSS_Added")), y = metadata.16s.1$shannon, fill = metadata.16s.1$Nitrogen)) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.75), dotsize = 0.65) +
scale_fill_manual(values=c('#b7a28c','#408080')) +
labs(title = 'Bacterial Evenness by Ecosystem', x = 'Ecosystem', y = 'Evenness', fill = 'Nitrogen') +
theme_classic() +
#theme(legend.position = "none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())
plot(bact.eco.nit.shannon)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

detach(metadata.16s.1)
###Fungi###
attach(metadata.ITS.1)
## The following objects are masked _by_ .GlobalEnv:
##
## Block, Collection_date, eco.nit, eco.prec, Nitrogen, Plot,
## Precipitation
Nitrogen<-factor(Nitrogen)
Precipitation <- as.factor(Precipitation)
Block <- as.factor(Block)
Plot <- as.factor(Plot)
Collection_date <- as.factor(Collection_date)
eco.prec <- as.factor(eco.prec)
eco.nit <- as.factor(eco.nit)
fung.aov.richness <- lmer(richness ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen + ( 1 | Block:Ecosystem), REML = FALSE, data = metadata.ITS.1)
## boundary (singular) fit: see ?isSingular
anova(fung.aov.richness, type = "III")
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Block 59.9 59.9 1 203 0.0634 0.801500
## Ecosystem 8863.5 8863.5 1 203 9.3822 0.002488 **
## Precipitation 3175.5 3175.5 1 203 3.3613 0.068209 .
## Nitrogen 89.9 89.9 1 203 0.0951 0.758053
## Collection_date 8574.7 1429.1 6 203 1.5127 0.175484
## Ecosystem:Precipitation 4453.3 4453.3 1 203 4.7139 0.031080 *
## Ecosystem:Nitrogen 1304.9 1304.9 1 203 1.3813 0.241258
## Ecosystem:Collection_date 11942.0 1990.3 6 203 2.1068 0.053957 .
## Precipitation:Collection_date 7921.5 1320.2 6 203 1.3975 0.217150
## Nitrogen:Collection_date 12097.0 2016.2 6 203 2.1341 0.050978 .
## Precipitation:Nitrogen 748.6 748.6 1 203 0.7924 0.374438
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fung.eco.prec.rich <- ggplot(data = fung.aov.richness, aes(x = factor(eco.prec, levels = c("Grass_Ambient","Grass_Reduced","CSS_Ambient", "CSS_Reduced")), y = metadata.ITS.1$richness, fill = metadata.ITS.1$Precipitation)) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center',
position=position_dodge(0.75), dotsize = 0.65) +
scale_fill_manual(values=c('#408080', '#FF8040')) +
labs(title = 'Fungal Richness by Ecosystem', x = 'Ecosystem', y = 'Richness', fill = 'Precipitation') +
theme_classic() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())
plot(fung.eco.prec.rich)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

fung.eco.prec.tuk.rich <- TukeyHSD(aov(formula = metadata.ITS.1$richness ~ as.character(metadata.ITS.1$eco.prec)))
fung.eco.prec.tuk.rich
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = metadata.ITS.1$richness ~ as.character(metadata.ITS.1$eco.prec))
##
## $`as.character(metadata.ITS.1$eco.prec)`
## diff lwr upr p adj
## CSS_Reduced-CSS_Ambient -18.877358 -36.41638 -1.3383342 0.0294106
## Grass_Ambient-CSS_Ambient -34.726415 -52.26544 -17.1873908 0.0000041
## Grass_Reduced-CSS_Ambient -36.071429 -53.94744 -18.1954220 0.0000026
## Grass_Ambient-CSS_Reduced -15.849057 -32.94844 1.2503233 0.0800862
## Grass_Reduced-CSS_Reduced -17.194070 -34.63893 0.2507867 0.0550132
## Grass_Reduced-Grass_Ambient -1.345013 -18.78987 16.0998433 0.9971695
fung.eco.nit.rich <- ggplot(data = fung.aov.richness, aes(x = factor(eco.nit, levels = c("Grass_Ambient","Grass_Added","CSS_Ambient", "CSS_Added")), y = metadata.ITS.1$richness, fill = metadata.ITS.1$Nitrogen)) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.75), dotsize = 0.65) +
scale_fill_manual(values=c('#b7a28c','#408080')) +
labs(title = 'Fungal Richness by Ecosystem', x = 'Ecosystem', y = 'Richness', fill = 'Nitrogen') +
theme_classic() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())
plot(fung.eco.nit.rich)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

fung.aov.shannon <- lmer(shannon ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen + ( 1 | Block:Ecosystem), REML = FALSE, data = metadata.ITS.1)
## boundary (singular) fit: see ?isSingular
anova(fung.aov.shannon, type = "III")
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Block 0.0026 0.00256 1 203 0.0152 0.902014
## Ecosystem 0.9103 0.91032 1 203 5.4137 0.020965 *
## Precipitation 0.2116 0.21159 1 203 1.2584 0.263286
## Nitrogen 0.0003 0.00033 1 203 0.0020 0.964757
## Collection_date 2.6120 0.43533 6 203 2.5889 0.019388 *
## Ecosystem:Precipitation 0.7041 0.70406 1 203 4.1871 0.042021 *
## Ecosystem:Nitrogen 0.3806 0.38064 1 203 2.2637 0.133992
## Ecosystem:Collection_date 3.2525 0.54208 6 203 3.2238 0.004786 **
## Precipitation:Collection_date 1.8596 0.30994 6 203 1.8432 0.092353 .
## Nitrogen:Collection_date 1.9503 0.32506 6 203 1.9331 0.077052 .
## Precipitation:Nitrogen 0.0168 0.01680 1 203 0.0999 0.752270
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fung.eco.prec.shannon <- ggplot(data = fung.aov.shannon, aes(x = factor(eco.prec, levels = c("Grass_Ambient","Grass_Reduced","CSS_Ambient", "CSS_Reduced")), y = metadata.ITS.1$shannon, fill = metadata.ITS.1$Precipitation)) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center',
position=position_dodge(0.75), dotsize = 0.65) +
scale_fill_manual(values=c('#408080', '#FF8040')) +
labs(title = 'Fungal Evenness by Ecosystem', x = 'Ecosystem', y = 'Evenness', fill = 'Precipitation') +
theme_classic() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())
plot(fung.eco.prec.shannon)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

fung.eco.prec.tuk.shannon <- TukeyHSD(aov(formula = metadata.ITS.1$shannon ~ as.character(metadata.ITS.1$eco.prec)))
fung.eco.prec.tuk.shannon
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = metadata.ITS.1$shannon ~ as.character(metadata.ITS.1$eco.prec))
##
## $`as.character(metadata.ITS.1$eco.prec)`
## diff lwr upr p adj
## CSS_Reduced-CSS_Ambient -0.20575542 -0.4467623 0.03525147 0.1235261
## Grass_Ambient-CSS_Ambient -0.41592122 -0.6569281 -0.17491433 0.0000766
## Grass_Reduced-CSS_Ambient -0.39190067 -0.6375381 -0.14626325 0.0003052
## Grass_Ambient-CSS_Reduced -0.21016580 -0.4451315 0.02479986 0.0975101
## Grass_Reduced-CSS_Reduced -0.18614525 -0.4258582 0.05356766 0.1870450
## Grass_Reduced-Grass_Ambient 0.02402055 -0.2156924 0.26373346 0.9938505
fung.eco.nit.shannon <- ggplot(data = fung.aov.shannon, aes(x = factor(eco.nit, levels = c("Grass_Ambient","Grass_Added","CSS_Ambient", "CSS_Added")), y = metadata.ITS.1$shannon, fill = metadata.ITS.1$Nitrogen)) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.75), dotsize = 0.65) +
scale_fill_manual(values=c('#b7a28c','#408080')) +
labs(title = 'Fungal Evenness by Ecosystem', x = 'Ecosystem', y = 'Evenness', fill = 'Nitrogen') +
theme_classic() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())
plot(fung.eco.nit.shannon)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

combined_rich_shan_bact_fun.precipitation <- ggarrange(bact.eco.prec.rich, bact.eco.prec.shannon, fung.eco.prec.rich, fung.eco.prec.shannon, labels = c("A.", "B.", "C.", "D."), ncol = 2, nrow = 2, widths = c(5,5)) # to remove legends and titles add "#" next to labs options of the individual plots.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
plot(combined_rich_shan_bact_fun.precipitation)

combined_rich_shan_bact_fun.nitrogen <- ggarrange(bact.eco.nit.rich, bact.eco.nit.shannon, fung.eco.nit.rich, fung.eco.nit.shannon, labels = c("A.", "B.", "C.", "D."), ncol = 2, nrow = 2, widths = c(5,5))
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
plot(combined_rich_shan_bact_fun.nitrogen)

detach(metadata.ITS.1)
###Bacteria###
median.avg.dist.16s <- avgdist(otu.clean.16s, sample = 1090, iterations = 1000, meanfun = median, dmethod = "bray")
## Warning in avgdist(otu.clean.16s, sample = 1090, iterations = 1000, meanfun
## = median, : The following sampling units were removed because they were below
## sampling depth: 14LRNDecember201716S, 14LRNSeptember201716S, 20RRNJune201716S,
## 20RRNSeptember201716S, 22LXXJune201716S, 22LXXSeptember201716S,
## 25LRXJune201716S, 32LXNDecember201716S, 35RRNJune201716S, 45LRXDecember201716S,
## 46RXXDecember201716S, 47LRNDecember201716S, 47RRXDecember201716S,
## 48LXXSeptember201716S, 48RXNDecember201716S, 4LXXAugust201616S, 4LXXJune201716S,
## 5LRNJune201716S, 8RXNJune201716S
bray.dist.16s <- as.data.frame(as.matrix(median.avg.dist.16s))
write.xlsx2(bray.dist.16s, file = "bray.dist.16s.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names = TRUE, append = FALSE)
NMDS.16s <- metaMDS(bray.dist.16s, distance = "bray", k = 3, noshare = FALSE, trymax = 400, trace = TRUE)
## Run 0 stress 0.1473617
## Run 1 stress 0.1470132
## ... New best solution
## ... Procrustes: rmse 0.01099409 max resid 0.109775
## Run 2 stress 0.1471217
## ... Procrustes: rmse 0.007527532 max resid 0.08744722
## Run 3 stress 0.147153
## ... Procrustes: rmse 0.009782592 max resid 0.1102538
## Run 4 stress 0.147169
## ... Procrustes: rmse 0.008878901 max resid 0.08772956
## Run 5 stress 0.1471381
## ... Procrustes: rmse 0.008596456 max resid 0.1099102
## Run 6 stress 0.1474069
## ... Procrustes: rmse 0.00947325 max resid 0.08680287
## Run 7 stress 0.1471773
## ... Procrustes: rmse 0.008628991 max resid 0.08740841
## Run 8 stress 0.1470301
## ... Procrustes: rmse 0.002391889 max resid 0.02444052
## Run 9 stress 0.1470599
## ... Procrustes: rmse 0.005266206 max resid 0.05515962
## Run 10 stress 0.1472593
## ... Procrustes: rmse 0.0115494 max resid 0.1112863
## Run 11 stress 0.1471472
## ... Procrustes: rmse 0.008735791 max resid 0.1100061
## Run 12 stress 0.147335
## ... Procrustes: rmse 0.01178672 max resid 0.109496
## Run 13 stress 0.1471293
## ... Procrustes: rmse 0.007341031 max resid 0.08740262
## Run 14 stress 0.1470149
## ... Procrustes: rmse 0.0001713161 max resid 0.001736122
## ... Similar to previous best
## Run 15 stress 0.1470132
## ... New best solution
## ... Procrustes: rmse 0.00212278 max resid 0.02365873
## Run 16 stress 0.147158
## ... Procrustes: rmse 0.008449116 max resid 0.1116341
## Run 17 stress 0.1472593
## ... Procrustes: rmse 0.01153399 max resid 0.1113063
## Run 18 stress 0.1472745
## ... Procrustes: rmse 0.01167785 max resid 0.1134128
## Run 19 stress 0.1472637
## ... Procrustes: rmse 0.01155497 max resid 0.111502
## Run 20 stress 0.146991
## ... New best solution
## ... Procrustes: rmse 0.001827964 max resid 0.0203683
## Run 21 stress 0.147257
## ... Procrustes: rmse 0.011816 max resid 0.1111026
## Run 22 stress 0.14701
## ... Procrustes: rmse 0.00176931 max resid 0.02012093
## Run 23 stress 0.147356
## ... Procrustes: rmse 0.009666349 max resid 0.08922747
## Run 24 stress 0.1471465
## ... Procrustes: rmse 0.006709326 max resid 0.0868903
## Run 25 stress 0.1472624
## ... Procrustes: rmse 0.01203293 max resid 0.1115047
## Run 26 stress 0.147162
## ... Procrustes: rmse 0.01013879 max resid 0.1127369
## Run 27 stress 0.1470135
## ... Procrustes: rmse 0.002892814 max resid 0.02483455
## Run 28 stress 0.1472321
## ... Procrustes: rmse 0.01145066 max resid 0.111831
## Run 29 stress 0.1469862
## ... New best solution
## ... Procrustes: rmse 0.002212889 max resid 0.02459004
## Run 30 stress 0.1472525
## ... Procrustes: rmse 0.01189783 max resid 0.1113853
## Run 31 stress 0.1472532
## ... Procrustes: rmse 0.01189261 max resid 0.1113517
## Run 32 stress 0.1471686
## ... Procrustes: rmse 0.008376918 max resid 0.08749068
## Run 33 stress 0.1470127
## ... Procrustes: rmse 0.002762964 max resid 0.02300611
## Run 34 stress 0.1471424
## ... Procrustes: rmse 0.009824468 max resid 0.1099773
## Run 35 stress 0.1470191
## ... Procrustes: rmse 0.003028171 max resid 0.02300531
## Run 36 stress 0.14707
## ... Procrustes: rmse 0.00526875 max resid 0.05637375
## Run 37 stress 0.1472604
## ... Procrustes: rmse 0.01187004 max resid 0.1114016
## Run 38 stress 0.1471355
## ... Procrustes: rmse 0.006614468 max resid 0.08691199
## Run 39 stress 0.1471404
## ... Procrustes: rmse 0.0090955 max resid 0.1102109
## Run 40 stress 0.1471245
## ... Procrustes: rmse 0.007358757 max resid 0.08761498
## Run 41 stress 0.1471608
## ... Procrustes: rmse 0.008310955 max resid 0.1103223
## Run 42 stress 0.1473535
## ... Procrustes: rmse 0.01084271 max resid 0.1099971
## Run 43 stress 0.1472248
## ... Procrustes: rmse 0.01112402 max resid 0.1112418
## Run 44 stress 0.1471302
## ... Procrustes: rmse 0.007576645 max resid 0.08783906
## Run 45 stress 0.1470569
## ... Procrustes: rmse 0.0046187 max resid 0.05439872
## Run 46 stress 0.1469986
## ... Procrustes: rmse 0.001149863 max resid 0.01225147
## Run 47 stress 0.1471684
## ... Procrustes: rmse 0.008967158 max resid 0.1072389
## Run 48 stress 0.1472329
## ... Procrustes: rmse 0.01143397 max resid 0.1111328
## Run 49 stress 0.1472266
## ... Procrustes: rmse 0.01131513 max resid 0.1112707
## Run 50 stress 0.1472576
## ... Procrustes: rmse 0.01204327 max resid 0.1115062
## Run 51 stress 0.1470195
## ... Procrustes: rmse 0.002260496 max resid 0.01956336
## Run 52 stress 0.1471483
## ... Procrustes: rmse 0.008977734 max resid 0.1106712
## Run 53 stress 0.1473286
## ... Procrustes: rmse 0.006317869 max resid 0.06006143
## Run 54 stress 0.1471239
## ... Procrustes: rmse 0.006817903 max resid 0.08728222
## Run 55 stress 0.147125
## ... Procrustes: rmse 0.007219689 max resid 0.08716617
## Run 56 stress 0.1470168
## ... Procrustes: rmse 0.002097874 max resid 0.02017475
## Run 57 stress 0.1471868
## ... Procrustes: rmse 0.00858156 max resid 0.08732222
## Run 58 stress 0.1472545
## ... Procrustes: rmse 0.01195251 max resid 0.1113325
## Run 59 stress 0.1471491
## ... Procrustes: rmse 0.009831946 max resid 0.1101909
## Run 60 stress 0.1471254
## ... Procrustes: rmse 0.006842557 max resid 0.08680285
## Run 61 stress 0.1472416
## ... Procrustes: rmse 0.008796869 max resid 0.08692925
## Run 62 stress 0.1470671
## ... Procrustes: rmse 0.005111794 max resid 0.05578186
## Run 63 stress 0.1471321
## ... Procrustes: rmse 0.00659307 max resid 0.08692297
## Run 64 stress 0.1473484
## ... Procrustes: rmse 0.01087259 max resid 0.1100298
## Run 65 stress 0.1472363
## ... Procrustes: rmse 0.01149527 max resid 0.1110132
## Run 66 stress 0.1471272
## ... Procrustes: rmse 0.007194208 max resid 0.08774139
## Run 67 stress 0.1472331
## ... Procrustes: rmse 0.01144492 max resid 0.111155
## Run 68 stress 0.1472372
## ... Procrustes: rmse 0.01099134 max resid 0.1114732
## Run 69 stress 0.1470017
## ... Procrustes: rmse 0.00258371 max resid 0.02365188
## Run 70 stress 0.1470111
## ... Procrustes: rmse 0.001886286 max resid 0.01936161
## Run 71 stress 0.1472254
## ... Procrustes: rmse 0.0112731 max resid 0.1112213
## Run 72 stress 0.1473473
## ... Procrustes: rmse 0.01099156 max resid 0.109975
## Run 73 stress 0.1472684
## ... Procrustes: rmse 0.01199933 max resid 0.1115151
## Run 74 stress 0.1471422
## ... Procrustes: rmse 0.006670685 max resid 0.0869103
## Run 75 stress 0.1473327
## ... Procrustes: rmse 0.01174141 max resid 0.1098759
## Run 76 stress 0.1472306
## ... Procrustes: rmse 0.01083578 max resid 0.1113879
## Run 77 stress 0.1471704
## ... Procrustes: rmse 0.008270334 max resid 0.08771603
## Run 78 stress 0.1473392
## ... Procrustes: rmse 0.01176678 max resid 0.1095753
## Run 79 stress 0.1471894
## ... Procrustes: rmse 0.007127218 max resid 0.09445938
## Run 80 stress 0.1473505
## ... Procrustes: rmse 0.01081933 max resid 0.1099293
## Run 81 stress 0.1470164
## ... Procrustes: rmse 0.002854988 max resid 0.02567395
## Run 82 stress 0.1471378
## ... Procrustes: rmse 0.007406443 max resid 0.0878771
## Run 83 stress 0.1473474
## ... Procrustes: rmse 0.0108633 max resid 0.1099446
## Run 84 stress 0.1472608
## ... Procrustes: rmse 0.01203316 max resid 0.1111735
## Run 85 stress 0.1470773
## ... Procrustes: rmse 0.004053615 max resid 0.04704992
## Run 86 stress 0.1474857
## ... Procrustes: rmse 0.008960734 max resid 0.08628311
## Run 87 stress 0.1470152
## ... Procrustes: rmse 0.002843881 max resid 0.02242059
## Run 88 stress 0.1469966
## ... Procrustes: rmse 0.002276933 max resid 0.02330802
## Run 89 stress 0.1472297
## ... Procrustes: rmse 0.01137699 max resid 0.1111843
## Run 90 stress 0.1471237
## ... Procrustes: rmse 0.006684047 max resid 0.08705333
## Run 91 stress 0.1471449
## ... Procrustes: rmse 0.009775783 max resid 0.1101108
## Run 92 stress 0.1472231
## ... Procrustes: rmse 0.01110888 max resid 0.1112124
## Run 93 stress 0.1472235
## ... Procrustes: rmse 0.01115106 max resid 0.1112283
## Run 94 stress 0.1470704
## ... Procrustes: rmse 0.003671529 max resid 0.04225642
## Run 95 stress 0.1471443
## ... Procrustes: rmse 0.01001536 max resid 0.1104978
## Run 96 stress 0.1473516
## ... Procrustes: rmse 0.01084375 max resid 0.1100857
## Run 97 stress 0.1471432
## ... Procrustes: rmse 0.009895986 max resid 0.1096512
## Run 98 stress 0.147237
## ... Procrustes: rmse 0.01144252 max resid 0.111118
## Run 99 stress 0.1470921
## ... Procrustes: rmse 0.005504368 max resid 0.05578434
## Run 100 stress 0.1471374
## ... Procrustes: rmse 0.006686398 max resid 0.08706045
## Run 101 stress 0.1471497
## ... Procrustes: rmse 0.01004118 max resid 0.1105268
## Run 102 stress 0.1471222
## ... Procrustes: rmse 0.006576593 max resid 0.08696867
## Run 103 stress 0.1473412
## ... Procrustes: rmse 0.01173908 max resid 0.1102461
## Run 104 stress 0.1471212
## ... Procrustes: rmse 0.00737984 max resid 0.08762399
## Run 105 stress 0.1471201
## ... Procrustes: rmse 0.007071948 max resid 0.08758696
## Run 106 stress 0.1472333
## ... Procrustes: rmse 0.01026069 max resid 0.1110523
## Run 107 stress 0.1473655
## ... Procrustes: rmse 0.01164181 max resid 0.1051452
## Run 108 stress 0.1472238
## ... Procrustes: rmse 0.01117147 max resid 0.1112131
## Run 109 stress 0.1472509
## ... Procrustes: rmse 0.01177104 max resid 0.1113409
## Run 110 stress 0.1471873
## ... Procrustes: rmse 0.008487483 max resid 0.08722651
## Run 111 stress 0.1473536
## ... Procrustes: rmse 0.01094308 max resid 0.1098447
## Run 112 stress 0.1472272
## ... Procrustes: rmse 0.01132252 max resid 0.1112118
## Run 113 stress 0.1469921
## ... Procrustes: rmse 0.002164316 max resid 0.02400088
## Run 114 stress 0.147131
## ... Procrustes: rmse 0.006668046 max resid 0.08717988
## Run 115 stress 0.1471793
## ... Procrustes: rmse 0.008121136 max resid 0.08704525
## Run 116 stress 0.1471438
## ... Procrustes: rmse 0.009907185 max resid 0.1100165
## Run 117 stress 0.1473549
## ... Procrustes: rmse 0.0108413 max resid 0.1099318
## Run 118 stress 0.1470648
## ... Procrustes: rmse 0.004396068 max resid 0.05314428
## Run 119 stress 0.1472324
## ... Procrustes: rmse 0.01143582 max resid 0.1111719
## Run 120 stress 0.1472248
## ... Procrustes: rmse 0.01123424 max resid 0.1112131
## Run 121 stress 0.1471215
## ... Procrustes: rmse 0.006601118 max resid 0.08702806
## Run 122 stress 0.1472351
## ... Procrustes: rmse 0.01135566 max resid 0.1111677
## Run 123 stress 0.1472358
## ... Procrustes: rmse 0.01154947 max resid 0.1111586
## Run 124 stress 0.14725
## ... Procrustes: rmse 0.01137833 max resid 0.1111902
## Run 125 stress 0.1473318
## ... Procrustes: rmse 0.01176556 max resid 0.1100493
## Run 126 stress 0.1469894
## ... Procrustes: rmse 0.0006347364 max resid 0.004565108
## ... Similar to previous best
## *** Solution reached
stressplot(NMDS.16s)

coordinates.16s <- data.frame(NMDS.16s$points[,1:3])
nmds.metadata.16s <- merge(coordinates.16s, metadata.16s.1, by.x = "row.names", by.y = "row.names")
nmds.metadata.16s <- data.frame(nmds.metadata.16s, col.names = TRUE, row.names = TRUE)
bacteria.nmds.plot <- ggplot(data = nmds.metadata.16s) +
aes(x = MDS1, y = MDS3, color = as.factor(Precipitation)) +
geom_point(aes(shape = as.factor(Ecosystem), size = as.factor(Ecosystem))) +
scale_color_manual(values = c('#408080', '#FF8040')) +
scale_shape_manual(values = c(17,16)) +
scale_size_manual(values = c(5,5)) +
geom_text(label = as.factor(nmds.metadata.16s$eco.prec), size = 1) +
theme_minimal() +
theme(legend.position = "none") +
labs(col = "Precipitation", shape = "Ecosystem") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text.x=element_blank(),
axis.ticks.x=element_blank(),axis.text.y=element_blank(),
axis.ticks.y=element_blank())
plot(bacteria.nmds.plot)

PERMANOVA.16s <- adonis(bray.dist.16s ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen, data = nmds.metadata.16s, strata = nmds.metadata.16s$Block, permutations = 999)
PERMANOVA.16s
##
## Call:
## adonis(formula = bray.dist.16s ~ Block + Ecosystem * Precipitation + Ecosystem * Nitrogen + Ecosystem * Collection_date + Collection_date * Precipitation + Collection_date * Nitrogen + Precipitation * Nitrogen, data = nmds.metadata.16s, permutations = 999, strata = nmds.metadata.16s$Block)
##
## Blocks: strata
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Block 1 3.206 3.2060 11.6324 0.04560 0.001 ***
## Ecosystem 1 4.019 4.0194 14.5836 0.05717 0.001 ***
## Precipitation 1 1.795 1.7948 6.5119 0.02553 0.001 ***
## Nitrogen 1 0.632 0.6316 2.2916 0.00898 0.001 ***
## Collection_date 6 7.456 1.2427 4.5090 0.10606 0.001 ***
## Ecosystem:Precipitation 1 1.283 1.2833 4.6563 0.01825 0.001 ***
## Ecosystem:Nitrogen 1 0.511 0.5109 1.8537 0.00727 0.003 **
## Ecosystem:Collection_date 6 2.805 0.4676 1.6965 0.03990 0.001 ***
## Precipitation:Collection_date 6 2.310 0.3851 1.3971 0.03286 0.001 ***
## Nitrogen:Collection_date 6 1.666 0.2777 1.0075 0.02370 0.337
## Precipitation:Nitrogen 1 0.524 0.5241 1.9015 0.00745 0.004 **
## Residuals 160 44.098 0.2756 0.62723
## Total 191 70.306 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###Fungi###
median.avg.dist.ITS <- avgdist(otu.clean.ITS, sample = 1064, iterations = 1000, meanfun = median, dmethod = "bray")
## Warning in avgdist(otu.clean.ITS, sample = 1064, iterations = 1000, meanfun
## = median, : The following sampling units were removed because they were
## below sampling depth: 20RRNJune2017ITS, 25LRXJune2017ITS, 27RXXAugust2016ITS,
## 27RXXSeptember2017ITS, 32LXNDecember2017ITS, 46RXXMarch2017ITS,
## 47RRXDecember2017ITS, 48LXXSeptember2017ITS, 5LRNDecember2017ITS,
## 8RXNMarch2017ITS
bray.dist.ITS <- as.data.frame(as.matrix(median.avg.dist.ITS))
write.xlsx2(bray.dist.ITS, file = "bray.dist.ITS.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names = TRUE, append = FALSE)
NMDS.ITS <- metaMDS(bray.dist.ITS, distance = "bray", k = 3, noshare = FALSE, trymax = 400, trace = TRUE)
## Run 0 stress 0.1034416
## Run 1 stress 0.102226
## ... New best solution
## ... Procrustes: rmse 0.02181735 max resid 0.2444773
## Run 2 stress 0.1034569
## Run 3 stress 0.1052886
## Run 4 stress 0.1022192
## ... New best solution
## ... Procrustes: rmse 0.0106754 max resid 0.08277034
## Run 5 stress 0.1066744
## Run 6 stress 0.1022795
## ... Procrustes: rmse 0.01290942 max resid 0.08527222
## Run 7 stress 0.1021602
## ... New best solution
## ... Procrustes: rmse 0.004342973 max resid 0.04197149
## Run 8 stress 0.1023053
## ... Procrustes: rmse 0.01315081 max resid 0.08762271
## Run 9 stress 0.1018855
## ... New best solution
## ... Procrustes: rmse 0.009392961 max resid 0.08142482
## Run 10 stress 0.1049468
## Run 11 stress 0.1024855
## Run 12 stress 0.1022005
## ... Procrustes: rmse 0.006606833 max resid 0.05844673
## Run 13 stress 0.1028988
## Run 14 stress 0.1034204
## Run 15 stress 0.1020033
## ... Procrustes: rmse 0.008730292 max resid 0.08205682
## Run 16 stress 0.1018634
## ... New best solution
## ... Procrustes: rmse 0.001273836 max resid 0.0101293
## Run 17 stress 0.1035167
## Run 18 stress 0.1036706
## Run 19 stress 0.1021602
## ... Procrustes: rmse 0.005574783 max resid 0.04095093
## Run 20 stress 0.1028066
## Run 21 stress 0.1047337
## Run 22 stress 0.101726
## ... New best solution
## ... Procrustes: rmse 0.007278107 max resid 0.08031383
## Run 23 stress 0.1075901
## Run 24 stress 0.1020344
## ... Procrustes: rmse 0.008923795 max resid 0.0820247
## Run 25 stress 0.103561
## Run 26 stress 0.1087161
## Run 27 stress 0.1094094
## Run 28 stress 0.1050336
## Run 29 stress 0.1036575
## Run 30 stress 0.1049014
## Run 31 stress 0.1033401
## Run 32 stress 0.1046354
## Run 33 stress 0.1026347
## Run 34 stress 0.1034715
## Run 35 stress 0.1069856
## Run 36 stress 0.1045293
## Run 37 stress 0.1032432
## Run 38 stress 0.1036727
## Run 39 stress 0.1017408
## ... Procrustes: rmse 0.002186305 max resid 0.01568619
## Run 40 stress 0.1022556
## Run 41 stress 0.1033054
## Run 42 stress 0.1047022
## Run 43 stress 0.1025933
## Run 44 stress 0.1057255
## Run 45 stress 0.1016655
## ... New best solution
## ... Procrustes: rmse 0.002711177 max resid 0.03583078
## Run 46 stress 0.1070481
## Run 47 stress 0.103016
## Run 48 stress 0.1025487
## Run 49 stress 0.1018629
## ... Procrustes: rmse 0.007001865 max resid 0.07879939
## Run 50 stress 0.1020274
## ... Procrustes: rmse 0.008085198 max resid 0.08109908
## Run 51 stress 0.1071755
## Run 52 stress 0.1036474
## Run 53 stress 0.1031518
## Run 54 stress 0.1022378
## Run 55 stress 0.1022273
## Run 56 stress 0.1025939
## Run 57 stress 0.1021349
## ... Procrustes: rmse 0.008153416 max resid 0.07793045
## Run 58 stress 0.1034201
## Run 59 stress 0.1040976
## Run 60 stress 0.1016764
## ... Procrustes: rmse 0.001591751 max resid 0.017169
## Run 61 stress 0.1036665
## Run 62 stress 0.1108792
## Run 63 stress 0.1030049
## Run 64 stress 0.1026753
## Run 65 stress 0.102592
## Run 66 stress 0.1020972
## ... Procrustes: rmse 0.006333918 max resid 0.06260569
## Run 67 stress 0.1019168
## ... Procrustes: rmse 0.007307465 max resid 0.08134799
## Run 68 stress 0.1128113
## Run 69 stress 0.1025593
## Run 70 stress 0.1054687
## Run 71 stress 0.109464
## Run 72 stress 0.1016627
## ... New best solution
## ... Procrustes: rmse 0.002280117 max resid 0.01837116
## Run 73 stress 0.102266
## Run 74 stress 0.1032318
## Run 75 stress 0.1019791
## ... Procrustes: rmse 0.008064308 max resid 0.08507382
## Run 76 stress 0.102034
## ... Procrustes: rmse 0.00815858 max resid 0.08194238
## Run 77 stress 0.1030287
## Run 78 stress 0.1045602
## Run 79 stress 0.1037836
## Run 80 stress 0.1052211
## Run 81 stress 0.1027585
## Run 82 stress 0.1029487
## Run 83 stress 0.1025479
## Run 84 stress 0.1025673
## Run 85 stress 0.1017324
## ... Procrustes: rmse 0.003368771 max resid 0.03510527
## Run 86 stress 0.1022964
## Run 87 stress 0.1036804
## Run 88 stress 0.108209
## Run 89 stress 0.1050465
## Run 90 stress 0.102855
## Run 91 stress 0.1020445
## ... Procrustes: rmse 0.007367707 max resid 0.06864915
## Run 92 stress 0.1024118
## Run 93 stress 0.1027368
## Run 94 stress 0.1055402
## Run 95 stress 0.1043254
## Run 96 stress 0.1052573
## Run 97 stress 0.1023046
## Run 98 stress 0.1057518
## Run 99 stress 0.1044599
## Run 100 stress 0.1049851
## Run 101 stress 0.1038992
## Run 102 stress 0.1025777
## Run 103 stress 0.1086398
## Run 104 stress 0.1052331
## Run 105 stress 0.1037097
## Run 106 stress 0.1017928
## ... Procrustes: rmse 0.003752577 max resid 0.0350888
## Run 107 stress 0.1053459
## Run 108 stress 0.1030664
## Run 109 stress 0.1026997
## Run 110 stress 0.1036137
## Run 111 stress 0.1019912
## ... Procrustes: rmse 0.004516958 max resid 0.04281677
## Run 112 stress 0.1023607
## Run 113 stress 0.1053249
## Run 114 stress 0.1018722
## ... Procrustes: rmse 0.007185418 max resid 0.0827704
## Run 115 stress 0.1049866
## Run 116 stress 0.1027422
## Run 117 stress 0.1016699
## ... Procrustes: rmse 0.002573967 max resid 0.01607748
## Run 118 stress 0.1032658
## Run 119 stress 0.1038003
## Run 120 stress 0.1037894
## Run 121 stress 0.1028438
## Run 122 stress 0.1022111
## Run 123 stress 0.107964
## Run 124 stress 0.1056293
## Run 125 stress 0.1020346
## ... Procrustes: rmse 0.007162788 max resid 0.06413547
## Run 126 stress 0.1033928
## Run 127 stress 0.1036479
## Run 128 stress 0.1045988
## Run 129 stress 0.1050291
## Run 130 stress 0.1046954
## Run 131 stress 0.1020148
## ... Procrustes: rmse 0.005192132 max resid 0.04981483
## Run 132 stress 0.1022318
## Run 133 stress 0.1021045
## ... Procrustes: rmse 0.00879418 max resid 0.08543031
## Run 134 stress 0.103443
## Run 135 stress 0.1033652
## Run 136 stress 0.104009
## Run 137 stress 0.1034414
## Run 138 stress 0.1026358
## Run 139 stress 0.1027934
## Run 140 stress 0.103136
## Run 141 stress 0.1045562
## Run 142 stress 0.104973
## Run 143 stress 0.1022534
## Run 144 stress 0.1032806
## Run 145 stress 0.10432
## Run 146 stress 0.1025979
## Run 147 stress 0.1025719
## Run 148 stress 0.1147042
## Run 149 stress 0.1027114
## Run 150 stress 0.1052497
## Run 151 stress 0.1026784
## Run 152 stress 0.1043505
## Run 153 stress 0.1025631
## Run 154 stress 0.1035082
## Run 155 stress 0.1047643
## Run 156 stress 0.1025771
## Run 157 stress 0.1020173
## ... Procrustes: rmse 0.007927123 max resid 0.08239197
## Run 158 stress 0.1026762
## Run 159 stress 0.105241
## Run 160 stress 0.1024549
## Run 161 stress 0.1019162
## ... Procrustes: rmse 0.005777645 max resid 0.06966694
## Run 162 stress 0.1022609
## Run 163 stress 0.1138861
## Run 164 stress 0.1041216
## Run 165 stress 0.1017235
## ... Procrustes: rmse 0.003172806 max resid 0.03560678
## Run 166 stress 0.1039789
## Run 167 stress 0.1025309
## Run 168 stress 0.1019256
## ... Procrustes: rmse 0.005060735 max resid 0.05253533
## Run 169 stress 0.1033446
## Run 170 stress 0.1139195
## Run 171 stress 0.1030149
## Run 172 stress 0.1018958
## ... Procrustes: rmse 0.005181505 max resid 0.05490719
## Run 173 stress 0.1051713
## Run 174 stress 0.1035176
## Run 175 stress 0.103605
## Run 176 stress 0.10508
## Run 177 stress 0.1030392
## Run 178 stress 0.1052827
## Run 179 stress 0.1027528
## Run 180 stress 0.1023023
## Run 181 stress 0.1018591
## ... Procrustes: rmse 0.004564508 max resid 0.05245591
## Run 182 stress 0.1017523
## ... Procrustes: rmse 0.004188954 max resid 0.03464671
## Run 183 stress 0.1072388
## Run 184 stress 0.1033439
## Run 185 stress 0.1025586
## Run 186 stress 0.1029229
## Run 187 stress 0.1025425
## Run 188 stress 0.1020394
## ... Procrustes: rmse 0.007627395 max resid 0.06618612
## Run 189 stress 0.1018584
## ... Procrustes: rmse 0.006902274 max resid 0.08236374
## Run 190 stress 0.1030269
## Run 191 stress 0.1029932
## Run 192 stress 0.10464
## Run 193 stress 0.1038516
## Run 194 stress 0.1018513
## ... Procrustes: rmse 0.006737391 max resid 0.08211791
## Run 195 stress 0.1025763
## Run 196 stress 0.1155397
## Run 197 stress 0.1053756
## Run 198 stress 0.1025485
## Run 199 stress 0.1033358
## Run 200 stress 0.105229
## Run 201 stress 0.1038723
## Run 202 stress 0.1019398
## ... Procrustes: rmse 0.005902901 max resid 0.05628087
## Run 203 stress 0.1034608
## Run 204 stress 0.1028819
## Run 205 stress 0.1053369
## Run 206 stress 0.1035027
## Run 207 stress 0.1051262
## Run 208 stress 0.1075342
## Run 209 stress 0.1026542
## Run 210 stress 0.1074899
## Run 211 stress 0.1045322
## Run 212 stress 0.1050699
## Run 213 stress 0.1032503
## Run 214 stress 0.1026424
## Run 215 stress 0.1037207
## Run 216 stress 0.1070578
## Run 217 stress 0.1023244
## Run 218 stress 0.1023371
## Run 219 stress 0.1081403
## Run 220 stress 0.1025895
## Run 221 stress 0.1038863
## Run 222 stress 0.1035582
## Run 223 stress 0.1061032
## Run 224 stress 0.1023903
## Run 225 stress 0.1027721
## Run 226 stress 0.1025444
## Run 227 stress 0.1051926
## Run 228 stress 0.1032486
## Run 229 stress 0.1025692
## Run 230 stress 0.1033395
## Run 231 stress 0.1027101
## Run 232 stress 0.1036359
## Run 233 stress 0.1033371
## Run 234 stress 0.1017805
## ... Procrustes: rmse 0.004208538 max resid 0.04196572
## Run 235 stress 0.1020692
## ... Procrustes: rmse 0.007853111 max resid 0.07078279
## Run 236 stress 0.1027522
## Run 237 stress 0.1022305
## Run 238 stress 0.1035704
## Run 239 stress 0.1073083
## Run 240 stress 0.1034984
## Run 241 stress 0.1052264
## Run 242 stress 0.1019462
## ... Procrustes: rmse 0.007695539 max resid 0.08465692
## Run 243 stress 0.1026462
## Run 244 stress 0.1034169
## Run 245 stress 0.1024427
## Run 246 stress 0.1057999
## Run 247 stress 0.103368
## Run 248 stress 0.1043746
## Run 249 stress 0.1025407
## Run 250 stress 0.1018606
## ... Procrustes: rmse 0.004734909 max resid 0.05428239
## Run 251 stress 0.1046786
## Run 252 stress 0.1024555
## Run 253 stress 0.1055564
## Run 254 stress 0.1019498
## ... Procrustes: rmse 0.005515581 max resid 0.05457991
## Run 255 stress 0.101737
## ... Procrustes: rmse 0.00339964 max resid 0.03563156
## Run 256 stress 0.1020452
## ... Procrustes: rmse 0.00844098 max resid 0.08612605
## Run 257 stress 0.1045068
## Run 258 stress 0.1017399
## ... Procrustes: rmse 0.004131214 max resid 0.02551933
## Run 259 stress 0.1022028
## Run 260 stress 0.101889
## ... Procrustes: rmse 0.00666658 max resid 0.07969948
## Run 261 stress 0.1043615
## Run 262 stress 0.1024206
## Run 263 stress 0.1052071
## Run 264 stress 0.1045263
## Run 265 stress 0.1020916
## ... Procrustes: rmse 0.008594322 max resid 0.07611796
## Run 266 stress 0.1049277
## Run 267 stress 0.1049517
## Run 268 stress 0.1054611
## Run 269 stress 0.1032181
## Run 270 stress 0.1047857
## Run 271 stress 0.1045819
## Run 272 stress 0.103215
## Run 273 stress 0.1034907
## Run 274 stress 0.105452
## Run 275 stress 0.1045778
## Run 276 stress 0.103836
## Run 277 stress 0.1054566
## Run 278 stress 0.1051983
## Run 279 stress 0.102321
## Run 280 stress 0.1035903
## Run 281 stress 0.1033324
## Run 282 stress 0.1016708
## ... Procrustes: rmse 0.002194432 max resid 0.02447812
## Run 283 stress 0.1045206
## Run 284 stress 0.1064861
## Run 285 stress 0.1037251
## Run 286 stress 0.1074868
## Run 287 stress 0.1027372
## Run 288 stress 0.1022857
## Run 289 stress 0.1050393
## Run 290 stress 0.1090944
## Run 291 stress 0.1037347
## Run 292 stress 0.1031956
## Run 293 stress 0.1049937
## Run 294 stress 0.103446
## Run 295 stress 0.1016554
## ... New best solution
## ... Procrustes: rmse 0.001915193 max resid 0.01460278
## Run 296 stress 0.1026395
## Run 297 stress 0.1032238
## Run 298 stress 0.108303
## Run 299 stress 0.105138
## Run 300 stress 0.1056164
## Run 301 stress 0.1045625
## Run 302 stress 0.1027069
## Run 303 stress 0.103958
## Run 304 stress 0.1025547
## Run 305 stress 0.1037276
## Run 306 stress 0.1017485
## ... Procrustes: rmse 0.003497165 max resid 0.02062046
## Run 307 stress 0.1019178
## ... Procrustes: rmse 0.004790139 max resid 0.05411592
## Run 308 stress 0.1020442
## ... Procrustes: rmse 0.005590636 max resid 0.05534946
## Run 309 stress 0.1060374
## Run 310 stress 0.1026752
## Run 311 stress 0.1018531
## ... Procrustes: rmse 0.007001605 max resid 0.08187253
## Run 312 stress 0.106804
## Run 313 stress 0.1044966
## Run 314 stress 0.1020695
## ... Procrustes: rmse 0.005723969 max resid 0.05691739
## Run 315 stress 0.1036345
## Run 316 stress 0.103325
## Run 317 stress 0.1019269
## ... Procrustes: rmse 0.004949187 max resid 0.05580597
## Run 318 stress 0.1022977
## Run 319 stress 0.1026207
## Run 320 stress 0.1025068
## Run 321 stress 0.1034426
## Run 322 stress 0.1037351
## Run 323 stress 0.1060372
## Run 324 stress 0.1020177
## ... Procrustes: rmse 0.007898818 max resid 0.08071515
## Run 325 stress 0.1035978
## Run 326 stress 0.1026315
## Run 327 stress 0.1022312
## Run 328 stress 0.1026271
## Run 329 stress 0.1028372
## Run 330 stress 0.1029809
## Run 331 stress 0.1071241
## Run 332 stress 0.1078917
## Run 333 stress 0.1025533
## Run 334 stress 0.104177
## Run 335 stress 0.1044862
## Run 336 stress 0.1050522
## Run 337 stress 0.1016668
## ... Procrustes: rmse 0.002204927 max resid 0.02284503
## Run 338 stress 0.1035263
## Run 339 stress 0.1050347
## Run 340 stress 0.105228
## Run 341 stress 0.1018703
## ... Procrustes: rmse 0.004501693 max resid 0.05519265
## Run 342 stress 0.1054364
## Run 343 stress 0.102541
## Run 344 stress 0.1043361
## Run 345 stress 0.1029396
## Run 346 stress 0.1049648
## Run 347 stress 0.1027477
## Run 348 stress 0.1034875
## Run 349 stress 0.1051458
## Run 350 stress 0.103004
## Run 351 stress 0.1040485
## Run 352 stress 0.1035455
## Run 353 stress 0.1022223
## Run 354 stress 0.1017381
## ... Procrustes: rmse 0.003034046 max resid 0.03617134
## Run 355 stress 0.1025491
## Run 356 stress 0.1108277
## Run 357 stress 0.1021067
## ... Procrustes: rmse 0.006369138 max resid 0.05852839
## Run 358 stress 0.1064937
## Run 359 stress 0.1088592
## Run 360 stress 0.1051294
## Run 361 stress 0.1018992
## ... Procrustes: rmse 0.004295914 max resid 0.05028761
## Run 362 stress 0.1022957
## Run 363 stress 0.1023214
## Run 364 stress 0.1062385
## Run 365 stress 0.1033258
## Run 366 stress 0.1027441
## Run 367 stress 0.1052934
## Run 368 stress 0.101762
## ... Procrustes: rmse 0.00351508 max resid 0.03640896
## Run 369 stress 0.1021703
## Run 370 stress 0.1030053
## Run 371 stress 0.103268
## Run 372 stress 0.1050429
## Run 373 stress 0.1035997
## Run 374 stress 0.103363
## Run 375 stress 0.1020397
## ... Procrustes: rmse 0.007650374 max resid 0.08010203
## Run 376 stress 0.1057458
## Run 377 stress 0.1031281
## Run 378 stress 0.1018803
## ... Procrustes: rmse 0.004253365 max resid 0.05252094
## Run 379 stress 0.1026492
## Run 380 stress 0.1049269
## Run 381 stress 0.102607
## Run 382 stress 0.1034051
## Run 383 stress 0.1022014
## Run 384 stress 0.104962
## Run 385 stress 0.1059308
## Run 386 stress 0.103671
## Run 387 stress 0.1021671
## Run 388 stress 0.1025342
## Run 389 stress 0.1036991
## Run 390 stress 0.1026817
## Run 391 stress 0.1037019
## Run 392 stress 0.101951
## ... Procrustes: rmse 0.005781329 max resid 0.05641397
## Run 393 stress 0.102302
## Run 394 stress 0.1099411
## Run 395 stress 0.1031918
## Run 396 stress 0.1034833
## Run 397 stress 0.1027881
## Run 398 stress 0.1054689
## Run 399 stress 0.1028868
## Run 400 stress 0.1059391
## *** No convergence -- monoMDS stopping criteria:
## 291: no. of iterations >= maxit
## 109: stress ratio > sratmax
stressplot(NMDS.ITS)

coordinates.ITS <- data.frame(NMDS.ITS$points[,1:3])
nmds.metadata.ITS <- merge(coordinates.ITS, metadata.ITS.1, by.x = "row.names", by.y = "row.names")
nmds.metadata.ITS <- data.frame(nmds.metadata.ITS, col.names = TRUE, row.names = TRUE)
fungi.nmds.plot <- ggplot(data = nmds.metadata.ITS) +
aes(x = MDS1, y = MDS3, color = as.factor(Precipitation)) +
geom_point(aes(shape = as.factor(Ecosystem), size = as.factor(Ecosystem))) +
scale_color_manual(values = c('#408080', '#FF8040')) +
scale_shape_manual(values = c(17,16)) +
scale_size_manual(values = c(5,5)) +
geom_text(label = as.factor(nmds.metadata.ITS$eco.prec), size = 1) +
theme_minimal() +
theme(legend.position = "none") +
labs(col = "Precipitation", shape = "Ecosystem") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text.x=element_blank(),
axis.ticks.x=element_blank(),axis.text.y=element_blank(),
axis.ticks.y=element_blank())
plot(fungi.nmds.plot)

PERMANOVA.ITS <- adonis(bray.dist.ITS ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen, data = nmds.metadata.ITS, strata = nmds.metadata.ITS$Block, permutations = 999)
PERMANOVA.ITS
##
## Call:
## adonis(formula = bray.dist.ITS ~ Block + Ecosystem * Precipitation + Ecosystem * Nitrogen + Ecosystem * Collection_date + Collection_date * Precipitation + Collection_date * Nitrogen + Precipitation * Nitrogen, data = nmds.metadata.ITS, permutations = 999, strata = nmds.metadata.ITS$Block)
##
## Blocks: strata
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Block 1 1.657 1.65720 7.2005 0.03203 0.001 ***
## Ecosystem 1 1.074 1.07426 4.6676 0.02076 0.001 ***
## Precipitation 1 0.615 0.61529 2.6734 0.01189 0.002 **
## Nitrogen 1 0.213 0.21280 0.9246 0.00411 0.495
## Collection_date 6 2.435 0.40575 1.7630 0.04705 0.001 ***
## Ecosystem:Precipitation 1 0.349 0.34851 1.5143 0.00674 0.070 .
## Ecosystem:Nitrogen 1 0.209 0.20875 0.9070 0.00403 0.502
## Ecosystem:Collection_date 6 2.225 0.37078 1.6110 0.04300 0.001 ***
## Precipitation:Collection_date 6 1.848 0.30796 1.3381 0.03571 0.022 *
## Nitrogen:Collection_date 6 1.404 0.23400 1.0167 0.02714 0.399
## Precipitation:Nitrogen 1 0.355 0.35459 1.5407 0.00685 0.075 .
## Residuals 171 39.356 0.23015 0.76067
## Total 202 51.738 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###Plant###
median.avg.dist.plant <- vegdist(otu.plant.2015.1, method = "bray")
bray.dist.plant <- as.data.frame(as.matrix(median.avg.dist.plant))
write.xlsx2(bray.dist.plant, file = "bray.dist.plant.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names = TRUE, append = FALSE)
NMDS.plant <- metaMDS(bray.dist.plant, distance = "bray", k = 3, noshare = FALSE, trymax = 400, trace = TRUE)
## Run 0 stress 0.1019617
## Run 1 stress 0.1022528
## ... Procrustes: rmse 0.03165542 max resid 0.08401119
## Run 2 stress 0.1019806
## ... Procrustes: rmse 0.01193097 max resid 0.03951897
## Run 3 stress 0.1003177
## ... New best solution
## ... Procrustes: rmse 0.03954919 max resid 0.141736
## Run 4 stress 0.1047205
## Run 5 stress 0.1045815
## Run 6 stress 0.1084601
## Run 7 stress 0.100317
## ... New best solution
## ... Procrustes: rmse 0.0005105433 max resid 0.001521973
## ... Similar to previous best
## Run 8 stress 0.1045886
## Run 9 stress 0.1019759
## Run 10 stress 0.1003173
## ... Procrustes: rmse 0.0003661874 max resid 0.001092883
## ... Similar to previous best
## Run 11 stress 0.1084616
## Run 12 stress 0.1003181
## ... Procrustes: rmse 0.0004177886 max resid 0.001676165
## ... Similar to previous best
## Run 13 stress 0.1052178
## Run 14 stress 0.1019863
## Run 15 stress 0.1019478
## Run 16 stress 0.1093034
## Run 17 stress 0.1084781
## Run 18 stress 0.1019414
## Run 19 stress 0.101944
## Run 20 stress 0.1003173
## ... Procrustes: rmse 0.0003675437 max resid 0.001147174
## ... Similar to previous best
## *** Solution reached
stressplot(NMDS.plant)

coordinates.plant <- data.frame(NMDS.plant$points[,1:3])
nmds.metadata.plant <- merge(coordinates.plant, metadata.plant.2015.1, by.x = "row.names", by.y = "row.names")
nmds.metadata.plant <- data.frame(nmds.metadata.plant, col.names = TRUE, row.names = TRUE)
plant.nmds.plot <- ggplot(data = nmds.metadata.plant) +
aes(x = MDS1, y = MDS2, color = as.factor(precipitation)) +
geom_point(aes(shape = as.factor(ecosystem), size = as.factor(ecosystem))) +
scale_color_manual(values = c('#408080', '#FF8040')) +
scale_shape_manual(values = c(17,16)) +
scale_size_manual(values = c(5,5)) +
geom_text(label = as.factor(nmds.metadata.plant$eco.prec), size = 1) +
theme_minimal() +
theme(legend.position = "right") +
labs(col = "Precipitation", shape = "Ecosystem") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text.x=element_blank(),
axis.ticks.x=element_blank(),axis.text.y=element_blank(),
axis.ticks.y=element_blank())
plot(plant.nmds.plot)

bray.dist.plant <- bray.dist.plant[(names(bray.dist.plant) %in% row.names(metadata.plant.2015.1)),]
bray.dist.plant <- bray.dist.plant[,(names(bray.dist.plant) %in% row.names(bray.dist.plant))]
matching.metadata.bray.plant <- metadata.plant.2015.1[c(match(sort(row.names(bray.dist.plant)), sort(row.names(metadata.plant.2015.1)))),]
PERMANOVA.plant <- adonis(bray.dist.plant ~ block + ecosystem*precipitation + precipitation*nitrogen + ecosystem*nitrogen, data = matching.metadata.bray.plant, strata = matching.metadata.bray.plant$block, permutations = 999)
PERMANOVA.plant
##
## Call:
## adonis(formula = bray.dist.plant ~ block + ecosystem * precipitation + precipitation * nitrogen + ecosystem * nitrogen, data = matching.metadata.bray.plant, permutations = 999, strata = matching.metadata.bray.plant$block)
##
## Blocks: strata
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## block 1 2.2688 2.26885 18.6187 0.26810 0.001 ***
## ecosystem 1 0.9055 0.90554 7.4311 0.10700 0.001 ***
## precipitation 1 1.4921 1.49207 12.2443 0.17631 0.001 ***
## nitrogen 1 0.0778 0.07777 0.6382 0.00919 0.659
## ecosystem:precipitation 1 0.6702 0.67015 5.4994 0.07919 0.006 **
## precipitation:nitrogen 1 0.0565 0.05651 0.4638 0.00668 0.795
## ecosystem:nitrogen 1 0.0672 0.06716 0.5511 0.00794 0.686
## Residuals 24 2.9246 0.12186 0.34559
## Total 31 8.4626 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#stacked bar plot of variance by factors
variance <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/variance_fig.xlsx")
variance <- as.data.frame(variance, col.names = TRUE, row.names = TRUE)
variance[1:3 , 1:3]
## factor community percent
## 1 Residuals Bacteria 59.86696
## 2 Residuals Fungi 81.81324
## 3 Residuals Plant 22.38277
c11 <- c("#cacccc",
"#4f4447",
"#697456",
"#91ab98",
"#504254",
"#408080",
"#54bfab",
"#306773",
"#809a97",
"#b48a6f",
"#9d7e72")
variance.plot <- ggplot(variance, aes(x = community, y = percent))+
geom_col(aes(fill = factor(factor, levels = c("Residuals", "Block", "Ecosystem", "Ecosystem_collection date","Collection_Date", "Precipitation", "Ecosystem_precipitation",
"Precipitation_collection date", "Precipitation_nitrogen", "Nitrogen", "Ecosystem_nitrogen"))), width = 0.7) +
theme_minimal() +
theme(axis.line.x=element_line(colour = "black"), axis.line.y=element_line(colour = "black")) +
scale_fill_manual(values = c11) +
labs(fill = "Factor (s)") +
#ggtitle("Variance Explained By Factor") + #Adds tittle and subtittle. Can modify p and r-squared values + title.
guides(size = "none") # makes the sizes guide dissapear
plot(variance.plot)

combined.nmds.variance.plots <- ggarrange(bacteria.nmds.plot, fungi.nmds.plot, plant.nmds.plot, variance.plot + font("x.text", size = 10), labels = c("A.", "B.", "C.", "D."), ncol = 2, nrow = 2, widths = c(5,5,5,4))
plot(combined.nmds.variance.plots)

ambient.water.input <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/ambient_water_input_figure.xlsx")
ambient.water.input <- as.data.frame(ambient.water.input, col.names = TRUE, row.names = TRUE)
ambient.water.input[1:2 , 1:3]
## month_year date ambient_water_input_mm
## 1 oct_2014 2014-10-01 0
## 2 oct_2014 2014-10-02 0
sapply(ambient.water.input, class)
## $month_year
## [1] "character"
##
## $date
## [1] "POSIXct" "POSIXt"
##
## $ambient_water_input_mm
## [1] "numeric"
ambient.water.input.1 <- ambient.water.input %>% mutate_if(is.POSIXt, as.Date)
ambient.water.input.2 <- ambient.water.input.1 %>% mutate_if(is.POSIXct, as.Date)
sapply(ambient.water.input.2, class)
## month_year date ambient_water_input_mm
## "character" "Date" "numeric"
reduced.water.input <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/reduced_water_input_figure.xlsx")
reduced.water.input <- as.data.frame(reduced.water.input, col.names = TRUE, row.names = TRUE)
reduced.water.input[1:2 , 1:3]
## month_year date reduced_water_input_mm
## 1 oct_2014 2014-10-01 0
## 2 oct_2014 2014-10-02 0
sapply(reduced.water.input, class)
## $month_year
## [1] "character"
##
## $date
## [1] "POSIXct" "POSIXt"
##
## $reduced_water_input_mm
## [1] "numeric"
reduced.water.input.1 <- reduced.water.input %>% mutate_if(is.POSIXt, as.Date)
reduced.water.input.2 <- reduced.water.input.1 %>% mutate_if(is.POSIXct, as.Date)
sapply(reduced.water.input.2, class)
## month_year date reduced_water_input_mm
## "character" "Date" "numeric"
pd = position_dodge(0.5)
water.fig <- ggplot() + geom_point(data = ambient.water.input.2, aes(x = date, y = ambient_water_input_mm), size=2, pch = 19, color = "#408080") +
geom_point(data = reduced.water.input.2, aes(x = date, y = reduced_water_input_mm), size=2, pch = 19, position = pd, color = "#a9805d") +
xlab('Dates') +
ylab('Precipitation (mm)') +
geom_line(data = ambient.water.input.2, aes(x = date, y = ambient_water_input_mm), color = "#408080") +
geom_line(data = reduced.water.input.2, aes(x = date, y = reduced_water_input_mm), position = pd, color = "#a9805d") +
scale_y_continuous(breaks = seq(0,70,by=5), expand = c(0, 0)) +
scale_x_date(breaks = date_breaks("month"), labels = date_format("%b"), expand = c(0, 0)) +
theme(panel.grid.major = element_line(colour = "black", linetype = "dotted"), panel.grid.minor = element_blank(), panel.grid.major.x= element_blank(),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.5)) +
scale_linetype_manual(values = c(1,1))
plot(water.fig)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
